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Abstract 


Efficient  feedback  control  algorithms  based  on  optimal  and  robust  control  theories 
have  been  formulated  and  tested  in  direct  numerical  simulations  of  turbulent  channel 
flow.  The  optimization  technique  used  is  based  solely  on  the  equations  governing  the 
fluid  flow  and  variations  of  a  mathematical  statement  of  the  control  objective,  without 
the  heuristic  procedures  normally  used  to  determine  effective  flow  control  algorithms. 
The  control  algorithms  tested  are  shown  to  be  extremely  effective,  and  a  host  of  new 
ideas  for  the  determination  of  simple,  implement  able,  effective  control  rules  for  turbulent 
flows  have  been  proposed  and  are  currently  still  under  investigation. 

The  problem  of  transition  control  via  optimal  and  robust  techniques  has  also  been 
studied  to  draw  parallels  between  the  linear  and  nonlinear  theories  on  problems  of  signif¬ 
icant  interest  in  fluid  mechanics.  Results  on  this  problem  have  also  been  quite  good  and 
clearly  demonstrate  how  the  control  theories  are  related.  With  this  insight,  an  important 
extension  of  the  concepts  of  robust  control  theory  to  nonlinear  problems  has  been  made. 


Accomplishments 

As  a  linear  “warm-up”  project,  optimal  and  robust  control  techniques  were  used  to 
effectively  control  small,  two-dimensional,  linearly  unstable  perturbations  to  a  laminar 
plane  channel  flow  at  Re  =  10,000.  The  outcome  was  control  rules  based  on  wall- 
information  only  which  were  highly  effective  at  stabilizing  the  flow  system,  and  is  dis¬ 
cussed  in  Part  A:  Optimal  and  robust  control  of  transition. 

The  application  of  control  theory  to  the  nonlinear  problem  of  turbulence  is,  of  course,  a 
much  greater  challenge.  The  model  problem  we  consider  is  turbulent  flow  in  a  plane  chan¬ 
nel  with  blowing  and  suction  distributed  over  the  walls  (as  an  idealization  of  boundary 
forcing  by  discrete  MEMS  actuators),  as  illustrated  in  figure  1.  In  the  case  of  nonlinear 
phenomena  such  as  turbulence,  iterative  approaches  must  be  used  based  on  local  lin¬ 
earizations  of  the  flow  state.  It  was  also  found  that  optimizations  must  be  performed 
over  finite  time  intervals  which  are  sufficiently  long  to  accurately  reflect  the  dynamical 
evolution  of  the  near- wall  flow.  After  some  difficulty,  an  extremely  effective  implementa¬ 
tion  of  the  optimal  control  technique  was  tested  which  reduced  the  drag  of  a  Re^  =  100 
channel  flow  by  approximately  50%.  This  far  exceeds  what  is  possible  using  heuristic 
techniques  in  the  same  flow.  The  calculations  also  illustrate  the  sensitivity  of  important 
integral  flow  quantities  (such  as  drag)  in  a  particular  flow  realization  to  small  modifica¬ 
tion  of  the  control  forcing;  thus,  a  valuable  new  tool  has  been  developed  which  niay  be 
used  to  identify  coherent  turbulent  structures  responsible  for  important  flow  characteris¬ 
tics  and,  more  importantly,  where  these  structures  are  most  sensitive  to  control  forcing. 
Results  are  discussed  in  Part  B:  Optimal  control  of  turbulence. 

Comparison  of  the  approaches  discussed  in  Parts  A  and  B  led  to  an  understanding 
of  how  the  concepts  from  linear  robust  control  theory  {i.e.  Hoo)  may  be  extended  to 
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(a)  Visualization  of  turbulent  channel  flow  at  Rcr  =  180.  Shaded  regions  indicate  coher¬ 
ent  structures  of  the  near-wall  turbulence.  Flow  is  from,  left  to  right,  walls  are  dark. 


(b)  Top  view  of  lower  wall,  which  is  covered  with  sensors  and  actuo.tors  in  one  possible 
configuration. 

Figure  1.  The  model  problem  studied  in  this  work  is  turbulent  flow  in  a  plane  channel.  Small 
amounts  of  blowing  and  suction  will  be  applied  through  the  computational  equivalent  of  closely 
spaced  holes  drilled  in  the  walls  in  response  to  these  turbulent  motions  in  a  manner  which 
reduces  drag.  This  is  a  (very  rough)  approximation  of  the  eventual  physical  implementation 
illustrated  in  (b). 

nonlinear  problems  in  a  consistent  manner.  System  robustness  is  achieved,  essentially, 
by  playing  the  “devils  advocate”  and  attempting  to  find  the  “best”  feedback  control 
in  the  presence  of  a  small  component  of  the  “worst-case”  disturbance  forcing  the  state 
equation.  This  type  of  forcing  leads  to  controls  which  are  less  prone  to  cause  instability  in 
the  system  in  the  presence  of  unmo defied  disturbances  at  the  price  of  a  slight  degradation 
of  performance  for  the  nominal  {i.e.  undisturbed)  plant.  Such  an  approach  is  easily  added 
to  the  optimal  control  algorithm  discussed  in  Part  B ,  and  is  put  in  a  rigorous  framework 
in  Part  C:  Robust  control  of  turbulence. 

The  final,  and  perhaps  most  important,  result  of  this  project  is  the  development  of  a 
technique  with  which  simple,  implementable  control  rules  may  be  rigorously  optimized 
with  similar  methods.  The  control  rules  under  consideration  in  this  portion  of  the  work 
are  based  on  fiow  information  which  may  be  obtained  with  flush  wall-mounted  sensors, 
and  determine  via  simple  feedback  rules  the  wall-normal  component  of  the  fluid  velocity 
distributed  on  the  walls.  The  techniques  of  optimal  and  robust  control  theory  are  used 
simply  to  optimize  the  unknown  coefficients  in  these  feedback  rules.  This  work  is  still 
under  investigation,  and  is  discussed  in  Part  D:  Optimization  of  practical  feedback 
rules  for  turbulence  control. 
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PART  A. 

Optimal  and  robust  control  of  transition 

Optimal  and  robust  control  theories  are  used  to  determine  feedback  control  rules  that 
effectively  stabilize  a  linearly  unstable  flow  in  a  plane  channel.  Wall  transpiration  (un¬ 
steady  blowing/suction)  with  zero  net  mass  flux  is  used  as  the  control.  Control  algo¬ 
rithms  are  considered  that  depend  both  on  full  flowfield  information  and  on  estimates  of 
that  flowfield  based  on  wall  skin-friction  measurements  only.  The  development  of  these 
control  algorithms  accounts  for  modeling  errors  and  measurement  noise  in  a  rigorous 
fashion;  these  disturbances  are  considered  in  both  a  structured  (Gaussian)  and  unstruc¬ 
tured  (“worst  case”)  sense.  The  performance  of  these  algorithms  is  analyzed  in  terms 
of  the  eigenmodes  of  the  resulting  controlled  systems,  and  the  sensitivity  of  individual 
eigenmodes  to  both  control  and  observation  is  quantified. 


1.  Introduction 

The  behavior  of  infinitesimal  perturbations  to  simple  laminar  flows  is  an  important 
and  well-understood  problem.  As  the  Reynolds  number  is  increased,  laminar  flows  often 
become  unstable  and  transition  to  turbulence  occurs.  The  effects  of  the  turbulence 
produced  in  such  flows  are  very  significant  and  often  undesirable,  resulting  in  increased 
drag  and  heat  transfer  at  the  flow  boundaries.  Thus,  a  natural  engineering  problem  is 
to  study  methods  of  flow  control  such  that  transition  to  turbulence  can  be  delayed  or 
eliminated. 

Transition  often  occurs  at  a  Reynolds  number  well  below  that  required  for  linear  in¬ 
stability  of  the  laminar  flow.  Orszag  &  Patera  (1983)  demonstrate  that  finite  amplitude 
two-dimensional  perturbations  can  highly  destabilize  infinitesimal  three-dimensional  per¬ 
turbations  in  the  flow.  Butler  &  Farrell  (1992)  show  that  the  non-orthogonality  of  the 
eigenmodes  of  subcriticai  flows  implies  that  perturbations  of  a  particular  initial  structure 
will  experience  large  amplification  of  energy  before  their  eventual  decay,  and  suggest  that 
such  amplification  can  sometimes  lead  to  flow  perturbations  large  enough  for  nonlinear 
instability  to  be  triggered.  Such  nonlinear  instabilities  can  lead  to  transition  well  below 
the  critical  Reynolds  number  at  which  linear  instability  occurs.  Results  such  as  these 
have  renewed  interest  in  the  control  of  the  small  (linear)  perturbations,  as  the  mitigation 
of  linear  perturbations  also  lessens  the  potency  of  these  nonlinear  “bypass”  mechanisms. 

A  firm  theoretical  basis  for  the  control  of  small  perturbations  in  viscous  shear  flows 
is  only  beginning  to  emerge.  An  important  step  in  this  direction  is  provided  by  Joslin 
et  ai  (1995)  and  Joshi,  Speyer,  &  Kim  (1996),  who  analyze  this  problem  in  a  closed- 
loop  framework,  in  which  the  dynamics  of  the  flow  system  together  with  the  controller 
are  examined.  Joslin  et  ai  (1995)  apply  optimal  control  theory  to  a  problem  related  to 
the  one  presented  here;  in  their  approach,  the  control  is  determined  through  an  adjoint 
formulation  requiring  full  flowfield  information.  Joshi,  Speyer,  &  Kim  (1996)  consider 
essentially  the  same  problem  analyzed  in  this  paper,  and  show  that  a  simple  constant 
gain  feedback  with  an  integral  compensator  may  be  used  in  a  single- input /single-output 
(SISO)  sense  to  stabilize  the  flow;  a  single  output  (the  appropriate  Fourier  component  of 
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the  streamwise  drag)  is  multiplied  by  some  scalar  K  and  summed  with  a  reference  signal 
to  determine  the  corresponding  component  of  the  control  velocity.  This  proportional 
approach  is  a  special  case  of  a  class  of  proportional-integral-derivative  (PID)  controllers, 
which  combine  terms  which  are  proportional,  integrals,  and  derivatives  of  a  scalar  output 
of  a  system. 

The  present  work  extends  these  analyses  to  rigorously  account  for  state  disturbances 
and  measurement  noise.  A  two-step  control  approach  is  used.  First,  a  state  estimate 
is  developed  from  a  (potentially  inaccurate)  model  of  the  flow  equations,  with  correc¬ 
tions  to  this  state  estimate  provided  by  (noisy)  flow  measurements  fed  back  through  an 
output  injection  matrix  L.  This  state  estimate  is  then  multiplied  by  a  feedback  matrix 
K  to  determine  the  control.  Potentially,  this  approach  can  yield  better  results  than  a 
PID  controller.  In  comparison  to  the  PID  approach,  the  present  approach  has  many 
more  parameters  in  the  control  law  (specifically,  the  elements  of  the  matrices  K  and  L), 
which  are  rigorously  optimized  for  a  clearly  defined  objective.  In  this  manner,  multiple- 
input /multiple-output  (MIMO)  systems  are  handled  naturally  and  the  controller  is  cou¬ 
pled  with  an  estimator  which  models  the  dynamics  of  the  system  itself. 

Though  a  PID  approach,  such  as  that  of  Joshi,  Speyer,  &  Kim  (1996),  is  sufficient  to 
stabilize  the  present  system,  it  is  the  authors’  judgement  that  application  of  modern 
control  theory  to  the  same  problem  is  a  timely  exercise.  Many  problems  in  fluid  me¬ 
chanics,  especially  those  involving  turbulence,  are  dominated  by  nonlinear  behaviour. 
In  such  problems,  the  linear  analysis  performed  in  this  paper  is  not  valid.  However, 
optimal  control  approaches,  which  use  full  state  information,  may  still  be  formulated 
(Abergel  &  Temam  1990)  and  performed  (Moin  &  Bewley  1995)  with  impressive  results. 
In  order  to  make  such  schemes  practical,  one  must  understand  how  to  account  for  distur¬ 
bances  in  a  rigorous  fashion  and  how  to  estimate  accurately  the  necessary  components 
of  the  state  (for  instance,  the  location  and  strength  of  the  near-wall  coherent  structures) 
based  on  limited  flow  measurements.  The  current  paper  makes  these  concepts  clear  in 
a  fluid-mechanical  sense,  albeit  for  a  linear  problem,  and  thus  provides  a  step  in  this 
development. 

The  controllers  and  estimators  used  in  this  work  are  determined  by  application  of  %2 
and  %oo  approaches.  These  techniques  have  recently  been  cast  in  a  very  compact  form 
by  Doyle  et  al  (1989),  and  are  well  suited  to  the  current  problem,  in  which  the  issue 
of  interest  is  the  ability  of  a  closed-loop  system  to  reject  disturbances  to  a  laminar  flow 
when  only  a  few  noisy  measurements  of  the  flow  are  available.  The  discussion  presented 
here  will  involve  some  tools  seen  often  in  the  controls  literature,  such  as  block  diagrams, 
which  are  not  in  common  use  by  the  fluids  community.  Such  tools  were  included  only 
after  careful  deliberation;  it  was  concluded  that  these  powerful  tools  are  essential  in 
making  this  development  clear,  and  are  thus  described  in  detail  when  used. 

In  §2,  we  derive  the  governing  equations  for  the  present  flow  stability  problem  and  cast 
these  equations  in  a  standard  notation,  which  makes  subsequent  application  of  control 
theory  straightforward.  In  §3,  the  control  problem  is  analyzed  in  terms  of  the  controlla¬ 
bility  and  observability  of  each  individual  eigenmode  of  the  system  developed  in  §2.  In 
§4,  the  control  approach  developed  in  Doyle  et  al.  (1989)  is  summarized  and  applied  to 
the  present  system.  In  this  control  approach,  two  Riccati  equations  describe  a  family  of 
7^2  ^-iid  PLoo  controllers,  which  take  into  account  structured  (Gaussian)  and  unstructured 
(“worst  case”)  disturbances.  Results  of  these  approaches  are  presented  in  §5,  and  §6 
presents  some  concluding  remarks. 
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PART  A.  Optimal  o.nd  robust  control  of  transition 

2.  Governing  equations 

This  chapter  derives  the  equations  governing  the  perturbations  to  a  laminar  channel 
flow  and  casts  them  in  a  form  to  which  standard  control  techniques  may  be  applied.  This 
familiar  discussion  is  presented  to  precisely  define  the  control  problem  under  considera¬ 
tion.  Readers  interested  only  in  how  the  control  techniques  are  applied  are  advised  to 
proceed  directly  to  §3. 


2.1.  Continuous  form,  of  flow  equations 

Consider  a  steady  plane  channel  flow  with  maximum  velocity  Uq  and  channel  half- width 
5.  Non-dimensionalizing  all  velocities  by  Uq  and  lengths  by  the  mean  velocity  profile  in 
the  streamwise  direction  (.t)  may  be  written  U{y)  =  l—y^  on  the  domain  y  E  [—1,1].  The 
equations  governing  small,  incompressible,  three-dimensional  perturbations  {u^v^w^p)  to 
the  mean  flow  U  are  given  by  linearized  Navier-Stokes  and  continuity 


d 

dx'^ 

t  dp 

AU^v  =-/  + 
ax 

<1 

(2.1a) 

d 

dx" 

It 

! 

+ 

i-A« 

Re 

(2.1b) 

d 

dx'^ 

dp 

‘—Aw 

Re 

(2.1c) 

du 

dv 

dw 

(2.2) 

dx 

where  A  =  jdx^  -h  fdy^  +  jdz^  is  the  Laplacian,  Re  =  UqS/v  is  the  Reynolds 
number,  u  is  the  kinematic  viscosity,  dot  (’)  denotes  9/5^,  and  prime  (')  denotes  d/dy. 
A  single  equation  for  the  normal  component  of  velocity  v,  found  by  taking  the  Laplacian 
of  (2.1b),  substituting  for  Ap  from  the  divergence  of  (2.1),  and  applying  (2.2),  is 

Av=  |-C/^A  +  ?7"£^  +  A(A/;?e)}  w.  (2.3a) 

The  equation  for  the  normal  component  of  vorticity  uj  —  du/dz  —  dw/dx^  found  by 
subtracting  dfdx  of  (2.1c)  from  djdz  of  (2.1a),  is 

61;=  { t;+  { -C/^  + A/iJe}  w.  (2.3b) 

As  the  domain  is  homogeneous  in  the  x  and  2:  directions,  we  may  Fourier  transform  the 
solution  such  that 


v{x,y,z,t)  =  ^  u(A;x,?y,^*;z,^)€xp[i(fca,.T-|-A:^z)] 
u){x,y,z,t)=  ^  u{kx,y,kz,t)exp[i{ka:X  A  kz  z)]. 


As  the  various  Fourier  modes  are  orthogonal  and  equations  (2.3a)  and  (2.3b)  are  linear, 
the  solution  for  each  wavenumber  pair  {k^^kz)  is  decoupled  and  obeys  the  equations 


Au  =  +  A(A/i?e)}  (2.4a) 

ci;  {-ikzU']  vA  {-ik^U  A /!:^lRe)  IV,  (2.4b) 


where  the  hat  accents  (")  have  been  dropped  for  notational  convenience  and  the  Laplacian 
now  takes  the  form  A  =  d'^fdy^  -  k^  -  k?^.  Equation  (2.4a)  is  the  (fourth  order)  Orr- 
Sommerfeld  equation  for  the  wall-normal  velocity  modes,  and  (2.4b)  is  the  (second  order) 
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equation  for  the  wall-normal  vorticity  modes.  Note  the  one-way  coupling  between  these 
two  equations.  Also  note  that,  from  any  solution  the  values  of  w,  and  p  may 

be  extracted  by  manipulation  of  the  above  equations  into  the  form 


(2.5a) 

(2.5b) 


Ajp  =  —2ikxU^v. 


(2.5c) 


Control  will  be  applied  at  the  wall  as  a  boundary  condition  on  the  wall-normal  component 
of  velocity  v.  The  boundary  conditions  on  u  and  w  are  no-slip  (ti  =  u;  =  0),  which  implies 
that  a;  =  0  and  (by  continuity)  (9u/5?/  =  0  on  the  wall. 

In  this  development,  it  is  assumed  that  an  array  of  sensors,  which  can  measure  stream- 
wise  and  spanwise  skin  friction,  and  actuators,  which  provide  wall-normal  blowing  and 
suction  with  zero  net  mass  flux,  are  mounted  on  the  walls  of  a  laminar  channel  flow. 
It  is  also  assumed  that  a  sufficient  number  of  sensors  and  actuators  are  installed  such 
that  individual  Fourier  components  of  wall  skin  friction  and  wall  transpiration  may  be 
approximated,  and  the  analysis  is  carried  through  for  a  particular  Fourier  mode. 


2.2.  Discrete  form  of  flow  equations 

The  continuous  problem  described  above  is  discretized  on  a  grid  of  iV  -h  1  Chebyshev- 
Gauss-Lobatto  points  such  that 

yi  —  cos(7r//iV)  for  0  ^  ^  iV. 

An  (AT-f- 1)  X  (A^  -{- 1)  matrix  V  may  be  expressed  (Canuto  et  al.  1988,  eqn.  2.4.31)  such 
that  the  derivative  of  u)  with  respect  to  y  on  the  discrete  set  of  AT  +  1  points  is  given  by 

iJ  =  ^  bj  and 

where  the  prime  (')  now  indicates  the  (partial)  derivative  of  the  discrete  quantity  with 
respect  to  y.  The  homogeneous  Neumann  boundary  condition  on  v  is  accomplished  by 
modifying  the  first  derivative  matrix  such  that 


4 


Differentiation  of  v  with  respect  to  y  is  then  given  by 

v”  =  ^  v”  ^  and  v””  = 

With  these  derivative  matrices,  it  is  straightforward  to  write  (2.4)  in  matrix  form. 
This  is  accomplished  by  first  expressing  the  matrix  form  of  (2.4)  on  all  A^-l- 1  collocation 
points  such  thatf 

V  =  Cv  (2.6a) 

cj  =  Cv-\-Sco,  (2.6b) 

where  £,  C,  and  S  are  (AT-f  1)  x  (N  1).  The  Dirichlet  boundary  conditions  are 
explicitly  prescribed  as  separate  “forcing”  terms.  To  accomplish  this,  decompose  £,  C, 

t  Note  that,  for  0?  the  matrix  form  of  the  LHS  of  (2.4a)  is  invertible,  so  the  form 

(2.6a)  is  easily  determined. 
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and  S  according  to 


where  An,  A21,  and  A22  are  (A^—  1)  x  {N  —  1)  and  ^n,  612,  621,  and  622  are  {N  —  1)  x  1. 
Noting  that  cjq  ==  (j^^n  =  0  by  the  no-slip  condition,  and  defining 


where  x  is  2{N  —  1)  x  1,  A  is  2{N  —  1)  x  2{N  —  1),  B  is  2(A^  —  1)  x  2,  and  it  is  2  x  1,  we 
may  express  (2.6)  in  the  standard  form 

X  —  Ax  -h  Bu.  (2.7) 

The  vector  x  is  referred  to  as  the  “state” ,  and  the  vector  u  is  referred  to  as  the  “control” . 
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where  ci ,  C2,  C3 ,  and  C4  are  1  x  {N 


Urti 


(  2/ml^ 
yrti2 

yms 

\2/m4  / 


1)  and  di,  ^2,  ds,  and  d4  are  1  x  1.  Finally,  defining 


a  Cl 

bcs 

fadi 

hdA 

ac2 

bc^ 

D  = 

ad2 

bd^ 

bci 

ac^ 

bdi 

ads 

hc2 

a  C4 

\bd2 

a  c?4  y 

where  is  4  X  1,  C  is  4  X  2{N  -  1),  and  D  is  4  x  2,  allows  us  to  express  ym  in  the 
standard  form  of  a  linear  combination  of  the  state  x  and  the  control  u 


Vm  =  Cx  +  Du. 

The  vector  y-m.  is  referred  to  as  the  “measurement” . 


(2.8) 


2.4.  Inner  products  and  norms 

For  this  paper,  the  inner  product  for  two  continuous  complex  functions  u  and  v  on  the 
domain  y  G  [— 1>  1]  is  defined  by 


(u,  V 


u*vC  dy, 


where  C(2/)  =  (1  “  2/^) 


and  the  star  (*)  denotes  the  complex  conjugate.  The  inner  product  for  functions  dis¬ 
cretized  on  the  collocation  points  yi  is  defined  by 


N 

m— 0 


where 


Cm=<^ 


TT 

m 

TT 

N 


m  =  0,  AT 


1  ^  m  ^  N  —  1. 


For  sufficiently  smooth  functions  u,v  on  a  sufficiently  large  number  N  of  Chebyshev- 
Gauss-Lobatto  grid  points  (Canuto  et  al.  1988),  the  discrete  inner  product  approximates 
the  continuous  inner  product,  (it,  v)n  ^  {u,  v)^.  Note  that,  for  two  discrete  vectors 
f  and  77  defined  only  on  the  interior  grid  points,  or  for  vectors  which  are  zero  at  the 
boundary  points  m  =  0,  AT,  the  inner  product  is  given  simply  by 


(^)  ‘n)N  - 


(2.9) 


where  star  (*)  applied  to  a  vector  denotes  conjugate  transpose.  The  norm  of  u,  denoted 
||u||,  is  defined  as  the  square  root  of  (u,u)  for  both  the  discrete  and  continuous  cases. 
Orthogonality  of  two  functions  implies  that  their  inner  product  is  zero. 


3.  Analysis  of  control  problem 

In  §2,  it  was  shown  that  the  equations  governing  small  perturbations  in  a  laminar 
channel  flow  may  be  expressed  in  the  standard  form 

x  =  Ax-{‘Bu  (3.1a) 

Vm  =  Cx  +  Du^  (3.1b) 

where  all  variables  are  complex  and  the  system  matrix  A  is  dense  and  non-self-adjoint. 
We  now  discuss  the  eigenmodes  of  A  and  identify  which  of  these  modes  may  be  modified 
by  the  control  u  and  which  may  be  detected  by  the  measurements  ym- 

It  has  been  shown  (Orszag  1971)  that,  for  i^e;$5772,  the  uncontrolled  problem  itself 
is  stable  and,  for  Re  >  5772,  weak  instability  is  seen  (though  most  of  the  eigenvalues 
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remain  stable),  with  the  greatest  instability  near  kx  =  1.0  and  =  0.0.  We  seek  a 
method  to  determine  the  control  u  which  stabilizes  the  system  in  a  manner  which  is 
robust  to  system  uncertainties.  To  simplify  our  discussion,  we  will  restrict  our  attention 
in  the  remainder  of  this  work  to  the  particular  case  Re  =  10,000,  A:®  =  !>  and  kz  =  0. 
Joshi,  Speyer,  &  Kim  (1996)  explore  the  {Re^  kx^  kz)  parameter  space  further. 

For  kz  =  0  (two-dimensional  perturbations),  C  =  0  in  (2.6),  entirely  decoupling  the 
(jj  eigenmodes  from  both  the  v  eigenmodes  and  from  the  control  u  =  (uq,  .  In  the 
language  of  control  theory,  the  lj  eigenmodes  are  thus  ‘‘uncontrollable”  by  the  control 
u,  (However,  it  is  also  seen  that  the  uj  eigenmodes  are  stable,  so  these  modes  will,  so  to 
speak,  “take  care  of  themselves”.)  Thus,  for  the  remainder  of  this  paper,  we  will  restrict 
our  attention  to  the  the  v  eigenmodes  according  to  system  (3.1)  with 


where  x  is  {N  —  1)  x  1,  A  is  {N  —  1)  x  (AT 


1),  .B  is  {N  —  1)  X  2,  and  '/x  is  2  x  1,  and 


^  ymi  \ 

( 

a  Cl 

\ 

fadi 

bd3\ 

Vm,  — 

2/m2 
ym3  1 

ac2 

6  Cl 

D  = 

ad2 

bdi 

bdi 

adz 

^  ymA  / 

\ 

bc2 

/ 

\bd2 

ad^J 

where  is  4  X  1,  (7  is  4  x  (AT  —  1),  and  B  is  4  x  2.  (All  the  constituent  matrices,  vectors, 
and  flow  measurements  are  described  in  the  previous  section.) 


3.1.  System  analysis 

We  now  address  whether  or  not  all  of  the  current  system’s  AT  —  1  eigenmodes  may  be 
controlled  by  the  m  =  2  control  variables,  and  whether  or  not  all  of  these  eigenmodes 
may  be  observed  with  the  p  =  4  measurements.  To  accomplish  this,  it  is  standard 
practice  to  consider  two  matrices  which  characterize  the  controllability  and  observability 
of  the  system  as  a  whole  (Lewis  1995).  These  are  the  system  controllability  Gramian  Lc 
of  (A,B)  and  the  system  observability  Gramian  Lo  of  (C,  A),  which  may  be  found  by 
solution  of 


ALc  +  LcA^  +  BB*  =0 
A*Lo-{-LoA  +  C^C  =  0. 

Note  that  stable  numerical  techniques  to  solve  equations  of  this  form,  referred  to  as 
Lyapunov  equations,  are  well  developed  (Kwakernaak  &  Sivan  1972). 

If  Lc  is  (nearly)  singular,  there  is  at  least  one  eigenmode  of  the  system  which  is  (nearly) 
unaffected  by  any  choice  of  control  w,  and  the  system  is  called  “uncontrollable”.  If  all 
uncontrollable  eigenmodes  are  stable,  and  a  controller  may  be  constructed  such  that  the 
dynamics  of  the  system  may  be  made  stable  by  the  application  of  control,  the  system  is 
called  “stabilizable” . 

Similarly,  if  Lq  is  (nearly)  singular,  there  is  at  least  one  eigenmode  of  the  system  which 
is  (nearly)  indiscernible  by  the  measurements  ym^  and  the  system  is  called  “unobserv¬ 
able”  .  If  all  unobservable  eigenmodes  are  stable,  and  an  estimator  may  be  constructed 
such  that  the  dynamics  of  the  error  of  the  estimate  may  be  made  stable  by  appropriate 
forcing  of  the  estimator  equation,  the  system  is  called  “detectable” . 

For  the  present  system,  the  smallest  eigenvalue  of  both  Lc  and  Lo  are  computed  to  be 
near  machine  zero,  indicating  that  the  present  system  as  derived  above  is  both  uncontrol- 
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lable  and  unobservable.  Gramian  analysis  can  not  identify  which  of  the  eigenmodes  are 
uncontrollable  or  unobservable,  however,  so  it  is  impossible  to  predict  from  this  analysis 
alone  whether  or  not  the  system  is  stabilizable  and  detectable.  For  this  reason,  we  now 
develop  a  method  to  determine  which  of  the  eigenmodes  of  a  system  may  be  affected  by 
the  control  u  and,  similarly,  which  eigenmodes  may  be  discerned  by  the  measurements 

2/m* 

3.2.  Individual  eigenmode  analysis 

We  will  now  make  use  of  the  modal  canonical  form  of  the  system  (3.1)  to  quantify  the 
sensitivity  of  each  eigenmode  of  A  to  both  control  and  observation  (Kailath  1980).  In 
order  to  clarify  the  derivation,  we  shall  examine  each  eigenmode  of  the  system  separately. 
Define  the  eigenvalues  and  the  right  and  left  eigenvectors,  and  77^,  of  A  such  that 


right  eigenvectors:  A^i  =  Xi^i 

left  eigenvectors:  77*  A  —  Xi  77^*, 

where  the  eigenvectors  are  normalized  such  that  ||^i||  =  ||77i||  =  1  for  all  z,  where  ||  *11  is 
defined  in  §2.4.  Assume  A  has  distinct  eigenvalues  (this  may  be  verified  for  the  present 
system  described  above).  Then  any  x  may  be  decomposed  as  a  linear  combination  of  the 
(independent  but  not  orthogonal)  right  eigenvectors  such  that 

x  =  '^ai^i.  (3.2a) 


DiflFerentiating  with  respect  to  time, 

x  =  '^diU.  (3.2b) 

i 

Also,  note  that  left  and  right  eigenvectors  corresponding  to  different  eigenvalues  are 
orthogonal 

(''?i)  Ci)  =  0  J  #  b  (3'3a) 

but  those  corresponding  to  the  same  eigenvalues  are  not 

(3.3b) 


3.2.1.  Definition  of  m,odal  control  sensitivity 
By  (3.1a)  and  (3.2),  we  have 

d  j  ^  ctj  gj + B  M 

i  i 

—  ^  ^ 
i 

Taking  the  inner  product  with  rjj  and  noting  (3.3a)  yields 

By  the  definition  of  the  inner  product  (2.9),  and  noting  (3.3b),  yields 


(3.4) 

(3.5) 


6i.j  —  Xj  ctj  -\- 


If  the  vector  B*77j  =  0,  then  Oj  =  Xjaj  for  any  u.  In  terms  of  equation  (3.2a),  the 
component  of  x  parallel  to  is  not  affected  by  the  control  w,  and  the  eigenmode  is  said 
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to  be  “uncontrollable” .  Further,  the  norm  of  the  coefficient  of  u 


(3.6) 


which  we  shall  call  the  control  sensitivity  of  mode  j,  is  a  quantitative  measure  of  the 
sensitivity  of  the  eigenmode  j  to  the  control  u.  Note  the  dependence  of  this  expression 
on  the  matrix  which  is  the  same  term  which  drives  the  Lyapunov  equation  for 

controllability  Gramian  Lc. 


3.2.2.  Definition  of  modal  observation  sensitivity 

By  (3.1b)  and  (3.2)  and  assuming,  for  the  moment,  that  w  =  0,  we  have 

Vm  —  ^  ^  C  . 

3 

If  the  vector  C  Cj  =  0,  then  ym,  will  not  be  a  function  of  aj.  In  terms  of  equation  (3.2a), 
the  component  of  x  parallel  to  does  not  contribute  to  the  measurements  and  the 
eigenmode  is  said  to  be  “unobservable”.  Further,  the  norm  of  C 

=  (3.7) 

which  we  shall  call  the  observation  sensitivity  of  mode  j,  is  a  quantitative  measure  of 
the  sensitivity  of  the  measurement  ym.  to  eigenmode  j.  Note  the  dependence  of  this 
expression  on  the  matrix  C*  C,  which  is  the  same  term  which  drives  the  Lyapunov 
equation  for  observability  Gramian  Lo- 

3.3.  Sensitivity  of  eigenm, odes  of  A  to  control  and  observation 

The  least  stable  eigenvalues  of  A  and  their  corresponding  control  and  observation  sensi¬ 
tivities  fj  and  Qj  are  tabulated  in  table  1.  Note  that  the  fourth  eigenmode  is  five  orders 
of  magnitude  less  sensitive  than  the  first  eigenmode  to  modifications  in  the  control.  In 
general,  those  modes  in  the  upper  branch  of  figure  la  (large  |^(A)|)  are  much  less  sensi¬ 
tive  to  control  than  those  in  the  lower  branch  (small  |^(A)|).  Near  the  intersection  of  the 
two  branches  (3i(A)  Ps:  —0.3),  the  control  sensitivity  is  maximum,  with  this  sensitivity 
decreasing  slowly  to  the  left  of  this  intersection  (K(A)  <  —0.3).  It  can  be  predicted  that 
the  eigenmodes  corresponding  to  the  largest  fj  may  be  affected  most  upon  application 
of  some  feedback  control  u. 

Note  that  the  flow  measurements  are  two  orders  of  magnitude  less  sensitive  to  the 
fourth  eigenmode  as  they  are  to  the  first  eigenmode.  It  can  be  predicted  that  the  state 
estimates  of  the  eigenmodes  corresponding  to  the  largest  gj  will  be  most  accurate  when 
estimating  the  state  based  on  the  measurements  in  the  presence  of  noise. 

An  important  observation  from  figure  lb  is  that  eigenvalues  in  the  upper  branch  of 
figure  la  have  corresponding  eigenvectors  with  variations  primarily  in  the  center  of  the 
channel,  and  are  thus  less  controllable  via  wall  transpiration  and  less  observable  via  wall 
measurements  than  eigenvalues  in  the  lower  branch.  This  observation  is  quantified  by 
reduced  values  of  fj  and  gj  for  these  modes  in  table  1. 

The  second  eigenvalue  computed,  at  A2  =  —0.0235  +  1.520  i  is  spurious.  Spurious 
eigenmodes  may  be  easily  identified  two  ways:  i)  the  eigenvalue  moves  significantly  when 
N  is  modified  slightly,  though  the  eigenvalues  reported  in  table  1  remain  converged,  and 
ii)  when  plotted,  spurious  modes  are  dominated  by  large  oscillations  from  grid  point 
to  grid  point  across  the  entire  domain,  though  converged  eigenmodes  are  well  resolved. 
Spurious  eigenmodes  are  expected  using  this  approach  and  may  be  disregarded. 
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ii 

i 

0.9| 

o.s; 

0.7; 

o.e! 

j 

0.5' 

i 

j 

0.4* 

0.3 

0.2 

0.1 


X  X 


-0.4  -0.3  -0.2  -0.1 


(a)  Least  stable  eigenvalues:  |^(Aj)|  versus  5R(Ay). 


(b)  Eigenvectors  corresponding  to  (left  to  right):  j  =  I  (unstable,  lower  branch), 
j  —  3  (stable,  upper  branch),  j  =  4  (stable,  upper  branch),  and  j  =  S  (stable, 
lower  branch),  plotted  as  a  function  of  y  from  the  lower  wall  (bottom)  to  the  upper 
wall  (top).  Real  component  of  eigenvector  is  shown  solid  and  imaginary  component 
dashed.  Corresponding  eigenvalues  are  reported  in  table  1. 


Figure  1.  Least  stable  eigenmodes  of  A  (no  control). 
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3 

Ai 

fj 

93 

1 

0.00373967-0.23752649?: 

0.266545 

102.61 

3 

-0.03516728  -  0.96463092  i 

0.000215 

72.85 

4 

-0.03518658  -  0.96464251  i 

0.000005 

1.45 

5 

-0.05089873  -  0.27720434?: 

0.026606 

347.98 

6 

-0.06320150  -  0.93631654?: 

0.000513 

81.39 

7 

-0.06325157-  0.93635178?: 

0.000021 

2.90 

8 

-0.09122274  -  0.90798305?: 

0.000931 

83.36 

9 

-0.09131286  -  0.90805633?: 

0.000056 

4.32 

10 

-0.11923285  -  0.87962729?: 

0.001587 

77.67 

11 

-0.11937073  -  0.87975570? 

0.000124 

5.37 

12 

-0.12450198  -  0.34910682?: 

0.171859 

69.50 

13 

-0.13822653  -  0.41635102?: 

0.037660 

252.09 

14 

-0.14723393  -0.85124584?: 

0.002833 

63.31 

15 

-0.14742560  -  0.85144938?: 

0.000268 

5.59 

16 

-0.17522868  -  0.82283504?: 

0.005581 

44.14 

38 

-0.32519719  -0.63610486?: 

5.659801 

0.78 

39 

-0.34373267  -  0.67764346  i 

4.685315 

0.64 

53 

-0.66286552  -  0.67027520?: 

0.259581 

11.58 

Table  1.  Least  stable  eigenmodes  of  A  (no  control)  and  the  associated  control  and  observation 
sensitivities.  Note  that  all  eigenvalues  agree  precisely  with  those  reported  by  Orszag  (1971). 
Calculation  used  Chebyshev  collocation  technique  with  N  —  140  in  quad  precision  (128  bits  per 
real  number).  The  second  eigenmode,  which  is  not  shown  here,  is  spurious  (see  text).  Note  that 
the  only  unstable  mode  {j  ===  1)  for  the  present  system  is  both  sensitive  to  the  control  u  and 
easily  detected  by  the  measurements  ijm- 


4.  Summary  of  7/2  and  Tioo  control  theories 

In  §2,  it  was  shown  that  the  equations  governing  small  perturbations  in  a  laminar 
channel  flow  may  be  expressed  in  the  standard  form 

X  —  Ax  4-  Bu  (4.1a) 

ym.  =  Cx  +  Du,  (4.1b) 

where  the  constituent  matrices  A,  B,  C,  and  D  were  summarized  and  discussed  in  §3. 
We  now  seek  a  simple  method  to  determine  a  control  u  based  on  the  measurements 
ym  to  force  the  state  x  towards  zero  in  a  manner  which  rigorously  accounts  for  state 
disturbances,  to  be  added  on  the  RHS  of  (4.1a),  and  measurement  noise,  to  be  added  on 
the  RHS  of  (4.1b).  Speciflcally,  we  will  consider  feedback  of  the  measurements  y^,  such 
that  a  state  estimate  x  is  first  determined  by  the  system  model 

X  =  Ax  +  Bu  —  u  (4.2a) 

fjm  =  Cx  +  Du,  (4.2b) 

u  =  Ti^ym  ~  i/m)>  (4.2c) 

then  this  state  estimate  is  used  to  produce  the  control 


u  =  }C{x), 


(4.3) 
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Equation  (4.1),  with  added  disturbance  terms  on  the  RHS,  is  referred  to  as  the  “plant”, 
(4.2)  is  referred  to  as  the  “estimator”,  and  (4.3)  is  referred  to  as  the  “controller”.  The 
problem  at  hand  is  to  compute  linear  time-invariant  (LTI)  functions  £  and  K  such  that  i) 
the  “output  injection”  term  u  forces  the  state  estimate  x  in  the  estimator  (4.2)  towards 
the  state  x  in  the  plant  (4.1),  and  ii)  the  control  u  computed  by  the  controller  (4.3)  forces 
the  state  x  towards  zero  in  the  plant  (4.1). 

The  flow  of  information  in  this  problem  is  illustrated  schematically  in  the  following 
block  diagram. 


disturbances 


measurements 

Vm 


i 


plant 


estimator 


state  estimate 

X 


controller 


control 

u 


The  plant,  which  is  forced  by  external  disturbances,  has  an  internal  state  x  which  cannot 
be  observed.  Instead,  a  few  noisy  measurements  ym  are  made,  and  with  these  measure¬ 
ments  an  estimate  of  the  state  x  is  determined.  This  state  estimate  is  then  fed  back  to 
through  the  controller  to  determine  the  control  u  to  apply  back  on  the  plant  in  order  to 
regulate  x  to  zero. 

We  will  now  demonstrate  how  to  apply  7^2  and  Tioo  control  theories  to  determine 
C  and  K.  (Note  that  we  will  redefine  several  variables  used  in  §2  to  derive  the  Orr- 
Sommerfeld  equation.  Considered  in  the  context  of  this  chapter,  this  should  present  no 
confusion.)  With  this  presentation,  one  set  of  control  equations,  involving  the  solution  of 
two  Riccati  equations,  describes  a  family  of  %2  and  Hoc  control  algorithms.  The  reader 
is  referred  to  Doyle  et  al.  (1989),  Dailey  et  al.  (1990),  and  Zhou,  Doyle,  &:  Glover  (1996) 
for  derivation  and  further  discussion  of  the  control  theories  summarized  here. 

4.1.  %2  control  theory 

4.1.1.  Optimal  control  (LQR) 

The  first  step  in  considering  the  system  (4.1)  is  to  consider  the  problem  with  no 
disturbances  and  measurements  which  identically  determine  full  information  about  the 
state,  so  that  x  =  x  {i.e.  no  estimation  of  the  state  is  necessary).  These  assumptions  are 
quite  an  idealization  and  can  rarely  be  accomplished  in  practice,  but  this  exercise  is  an 
important  step  to  determine  the  best  possible  system  performance.  It  is  for  this  reason 
that  the  controller  in  this  limit  is  referred  to  as  optimal.  Under  these  assumptions 
about  the  system,  the  objective  of  the  optimal  controller,  of  the  form  in  (4.3),  is  to 
regulate  {i.e.  return  to  zero)  some  measure  of  the  flow  perturbation  x  from  an  arbitrary 
initial  condition  as  quickly  as  possible  without  using  excessive  amounts  of  control  forcing. 
Mathematically,  a  cost  function  for  this  problem  may  thus  be  expressed  as 


(4.4) 
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The  term  involving  ||x'|p  is  a  measure  of  the  flow  perturbation  x  integrated  over  the 
time  period  over  which  this  perturbation  decays,  which  is  taken  as  f  G  [0,oo).  The 
term  involving  is  an  expression  of  the  magnitude  of  the  control.  These  two  terms  are 
weighted  together  with  a  scalar  ,  which  represents  the  price  of  the  control.  This  quantity 
is  small  if  the  control  is  ‘^cheap”  (which  generally  results  in  larger  control  magnitudes), 
and  large  if  applying  the  control  is  “expensive” .  As  the  state  equation  is  linear,  the  cost 
quadratic,  and  the  control  objective  regulation,  this  controller  is  also  referred  to  as  a 
linear  quadratic  regulator  (LQR). 

The  mathematical  statement  of  the  present  control  problem,  then,  is  the  minimization 
of  Jlqr-  This  results  in  regulation  of  x  without  excessive  use  of  control  effort.  Note 
that  minimization  of  Jlqb.  is  equivalent  to  minimization  of  the  integral  of  where 


2;  = 


and  where  Q  is  a  diagonal  matrix  with  diagonal  entries  Qjj  =  tt/AT,  as  required  by  the 
deflnition  of  the  norm  in  §2.4.  In  order  to  arrive  at  a  form  which  is  easily  generalized  in 
later  sections,  define 


B2  =  B 


Di2  =1 


For  notational  convenience,  the  state  equation  (4.1a)  will  be  considered  as  “forced”  with 
a  right  hand  side  forcing  term  r  which  shall  be  set  to  zero,  as  this  regulation  problem 
simply  drives  the  state  towards  zero  without  external  command  input.  The  state  equation 
(4.1a),  the  performance  measure  2;,  and  the  state  estimate  x  then  may  be  written 


X  =  Ax  +  r  B2  u 

(4.5a) 

z  =  Cix  +  D12U 

(4.5b) 

X  =  X. 

(4.5c) 

The  optimal  controller  Klqb.  is  sought  to  relate  the  control  u  to  the  (precise)  state 
estimate  x.  Control  is  applied  to  modify  the  evolution  of  the  state  x  such  that  the  cost 
Jlqr{z)  is  minimized.  The  important  matrices  of  the  system  described  by  (4.5)  may  be 
summarized  in  the  shorthand  form 


X 

r 

u 

X 

A 

I 

B2 

riQB  =  Z 

Cl 

0 

Di2 

X 

I 

0 

0 

The  flow  of  information  is  represented  by  the  block  diagram 


Vlqr 


r  =  0 


K^lqr 


where  Vlqr  is  the  flow  system  given  by  (4.5)  and  JClqr  is  the  optimal  controller,  which  is 
still  to  be  determined.  Note  that  the  command  input  is  r  =  0  and  there  are  no  disturbance 
inputs;  the  task  of  the  control  u  is  simply  to  regulate  the  state  x  from  nonzero  initial 
conditions  back  to  zero.  The  state  x  =  x  is  fed  back  through  the  controller  JClqr  to 
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control  the  system.  The  system  output  may  be  used  to  monitor  the  performance  of 
the  system. 

Given  this  general  setup,  a  Hamiltonian  is  defined  such  that 


H2 


(  A 

\-ClC,  -A*  )■ 


(4.6a) 


As  shown  in  Doyle  et  al.  (1989),  the  Hermetian  positive-definite  solution  X2  to  the  alge¬ 
braic  Riccati  equation  defined  by  this  Hamiltonian 


A*  X2  +  X2  A  -  X2  {B2  b; )  X2  +  (C*  Cl )  =  0,  (4.6b) 


denoted  X2  =  Ric(iJ2),  then  yields  the  optimal  LTI  state  feedback  matrix 

K2  =  -b;x2. 


(4.6c) 


The  optimal  LTI  controller  Klqr  then  given  simply  by 

(4.7) 

This  controller  minimizes  z*zdt  in  a  system  with  no  disturbances  and  arbitrary  initial 
conditions.  Note  that  standard  numerical  techniques  to  solve  equations  of  the  form  (4.6b) 
are  well  developed  (Laub  1991). 


u  =  K2  X 


4.1.2.  Kalman-Bucy  filter  (KBF) 

When  there  are  disturbances  to  the  system,  and  thus  the  state  is  not  precisely  known, 
the  state  (or  some  portion  thereof)  must  first  be  estimated,  then  the  control  determined 
based  on  this  state  estimate.  The  Kalman-Bucy  filter,  of  the  form  (4.2),  accomplishes  the 
required  state  estimation  by  assuming  that  the  state  disturbances  and  the  measurement 
noise  are  uncorrelated  white  Gaussian  processes.  To  accomplish  this,  we  introduce  two 
zero-mean  white  Gaussian  processes  wi  and  W2  with  covariance  matrices  E[wlwi]  =  /, 
E[w2W2]  =  /,  where  E[‘]  denotes  the  expectation  value.  With  these  new  disturbance 
signals,  and  with  Gi  defined  as  the  square  root  of  the  covariance  of  the  disturbances  to 
the  state  equation  and  G2  defined  as  the  square  root  of  the  covariance  of  measurement 
noise,  the  system  (4.1)  takes  the  form 

X  =  Ax  +  GiWi  +  Bu  (4.8a) 

ym  —  Gx  -f  G2W2  +  Du.  (4.8b) 

The  objective  of  the  Kalman-Bucy  filter  is  to  estimate  the  state  x  as  accurately  as  possible 
based  solely  on  the  measurements  ym.>  Put  another  way,  the  Kalman-Bucy  filter  attempts 
to  regulate  the  norm  of  the  state  estimation  error  xe  to  zero,  where 

X  E  =  X  —  X 

and  where  the  state  estimate  x  shall  be  determined  by  a  filter  of  the  form  (4.2).  Mathe¬ 
matically,  a  cost  function  for  this  problem  may  thus  be  expressed  as 

JKBF  =  E[\\zEfl  (4.9) 

where  ze  ^  xe  for  notational  convenience.  (As  Gaussian  disturbances  'll?!  and  W2  con¬ 
tinually  drive  this  system,  an  integral  on  ^  €  [0,oo),  as  used  to  define  Jlqr^  is  not 
convergent  for  this  problem,  and  the  expectation  value  is  the  relevant  measure.) 

The  mathematical  statement  of  the  present  control  problem,  then,  is  the  minimization 
of  JkbF'  This  results  in  a  “best  possible”  estimate  of  the  state  x.  In  order  to  arrive  at 
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a  form  which  is  easily  generalized  in  later  sections,  assume  G2  is  nonsingular  and  define 
5i  =  (Gi  0)  C2  =  GpC  D21  =  (0  I) 

and  the  vector  of  disturbances 

Also,  define  new  “observation”  vectors  y  and  y  by  a  simple  change  of  variables  such  that 

y  =  Gp{ym-  Du)  y  =  Gp (ym  -Du). 

Note  that  this  change  of  variables  does  not  represent  any  real  limitation,  for  whenever  any 
flow  measurement  ym.  is  made  in  a  physical  implementation,  the  control  u  at  that  moment 
is  also  known,  so  the  observation  y  is  easily  determined  from  the  flow  measurement  ym.- 
With  this  change  of  variables,  (4.8b)  and  (4.2b)  may  be  expressed  as 

y  =  C2X  D2\  w  (4.10a) 

y  =  C2X.  (4.10b) 

As  we  are  developing  the  equations  for  an  estimator,  it  is  appropriate  now  to  examine  the 
equations  for  the  state  estimation  error  xe  and  the  ouput  estimation  error  yE  =  y  —  y- 
Subtracting  (4.2a)  from  (4.8a)  and  (4.10b)  from  (4.10a)  yields  the  system 


xe  —  Axe  -\-Biw-{-u  (4.11a) 

ze  =  XE  (4.11b) 

ys  =  C2XE  +  D21W.  (4.11c) 

The  Kalman-Bucy  filter  Ckbf  is  sought  to  relate  the  output  injection  term  u  to  the 
output  estimation  error  yE-  The  extra  term  u  is  applied  in  these  model  equations  to 
control  the  evolution  of  the  state  estimation  error  xe  such  that  the  cost  Jkbf{^e)  is 
minimized  in  the  presence  of  Gaussian  disturbances  to.  The  important  matrices  of  the 
system  described  by  (4.11)  may  be  summarized  in  the  shorthand  form 


XE 

w 

u 

XE 

r  A 

I 

Vkbf  = 

I 

0 

0 

VE 

.  C2 

D21 

0 

The  flow  of  information  is  represented  by  the  block  diagram 

Ze  I - 1  w 


VE 


Vkbf 


T'Kbf 


where  Vkbf  is  the  flow  system  given  by  (4.11)  and  Ckbf  is  the  Kalman-Bucy  filter,  which 
is  still  to  be  determined.  This  system  accounts  for  Gaussian  disturbances  w  and  noisy 
observations  yE  of  the  system,  which  are  fed  back  through  the  filter  Ckbf  to  produce 
the  state  estimate.  The  system  output  ze  may  be  used  to  monitor  the  performance  of 
the  system.  Note  the  striking  similarity  of  the  structure  of  Vkbf  to  the  structure  of 
the  conjugate  transpose  of  Vlqr-  For  this  reason,  these  two  problems  are  referred  to  as 
“duals” ,  and  their  solutions  are  closely  related. 
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Given  this  general  setup,  another  Hamiltonian  is  defined  such  that 


J2  = 


-Bi  B* 


-c;c2\ 

-A  ) 


(4.12a) 


As  shown  in  Doyle  et  al  (1989),  the  Hermetian  positive-definite  solution  Y2  to  the  alge¬ 
braic  Riccati  equation  defined  by  this  Hamiltonian 


A  ^2  +  1^2  A*  -  Y2  (C*  C2)  Y2  +  (Hi  )  =  0, 

denoted  Y2  =  Ric(J2),  then  yields  the  LTI  estimator  feedback  matrix 

L2  =  -Y2c;, 

The  LTI  Kalman-Bucy  filter  Ckbf  is  then  simply  given  by 

u  =  L2  VE, 


(4.12b) 


(4.12c) 


and  thus  the  complete  state  estimator  is  given  by 


X  =  Ax  -\r  B2U  ~  1/2  {y  —  C2  x) 


(4.13) 


This  estimator  minimizes  E[\\x  —  .x*|p]  in  a  system  with  Gaussian  disturbances  in  the 
state  equation  and  Gaussian  noise  in  the  measurements. 


4.1.3.  %  control  (LQG  =  LQR  -h  KBF) 

An  estimator/controller  of  the  form  (4.2) — (4.3)  for  the  complete  system  described  by 
(4.8)  with  Gaussian  disturbances  may  now  be  constructed.  The  objective  of  the  control 
is  to  minimize 


J2^E[\\xf  +fu*u\,  (4.14) 

where  ||  -  ||  denotes  the  “2-norm”  as  defined  in  §2.5.  Note  that  minimization  of  J2  is 
equivalent  to  minimization  of  the  expectation  value  of  where 

«  j’ 

and  Q  is  a  diagonal  matrix  with  diagonal  entries  Qjj  =  tt/AT  as  required  by  the  definition 
of  the  norm  in  §2.4.  As  the  control  objective  is  the  minimization  of  the  expectation  value 
of  the  square  of  a  2-norm,  this  type  of  estimator /controller  is  referred  to  as  7i2.  As  the 
state  equation  is  linear,  the  cost  quadratic,  and  the  disturbances  Gaussian,  this  type  of 
estimator /controller  is  also  referred  to  as  linear  quadratic  Gaussian  (LQG). 

Combining  the  notation  developed  in  the  previous  two  sections 


Bi  =  {Gi  0)  Cl  =  (*^'0  Di2  = 

B2  =  B  C2^G;^C  £>2i  =  (0  /), 

with  the  vector  of  disturbances  w  and  the  observation  vectors  y  and  y  defined  such  that 

y  =  G2  {ym  —  Du) 


_ 


V  =  G2  ^{Vm  -Du), 
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the  system  (4.8)  and  the  control  objective  for  the  minimization  of  J2  take  the  form 

x  —  Ax  +Biw-\-B2U  (4.15a) 

2;  =  Cix-\-  D12U.  (4.15b) 

y  =  C2X  +  D21W.  (4.15c) 


An  ^^2  estimator/controller  is  sought  to  relate  the  observations  y  to  the  control  which 
is  applied  to  control  the  evolution  of  the  state  x  such  that  the  cost  J2{z)  is  minimized 
in  the  presence  of  Gaussian  disturbances  w.  The  important  matrices  of  the  complete 
system  described  by  (4.15)  may  be  summarized  in  the  shorthand  form 


X 

w 

u 

X 

A 

Bi 

B2 

P  complete  —  Z 

Cl 

0 

Di2 

y 

.  C2 

D21 

0 

Similarly,  the  important  matrices  of  the  model  plant  in  an  estimator  of  the  form  (4.2) 
may  be  represented  in  the  shorthand  form 


X  u  u 


X 

A 

-I 

B2 

P model  ~ 

I 

0 

0 

V 

.  C2 

0 

0 

With  these  representations,  the  flow  of  information  is  represented  by  the  block  diagram 


system  output  disturbances 


(controller) 


The  plant,  which  is  forced  by  external  disturbances  w;,  has  an  internal  state  x  which 
cannot  be  observed.  Instead,  a  few  noisy  observations  y  are  made,  and  with  these 
observations  an  estimate  of  the  state  x  is  determined.  This  state  estimate  is  then  fed 
back  to  through  the  controller  to  determine  the  control  u  to  apply  back  on  the  plant  in 
order  to  regulate  x  to  zero.  The  system  output  2:  may  be  used  to  monitor  the  performance 
of  the  system. 

The  remarkable  result  from  control  theory  (Lewis  1995)  is  that  the  I-L2  estimator/controller 
of  the  form  illustrated  in  the  above  block  diagram  which  minimizes  J2  for  this  system  is 
formed  by  simple  combination  of  the  optimal  controller  and  the  Kalman-Bucy  filter  such 
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X  =  Ax  B2U-  L2  {y  -  C2  x) 
u  =  K2  X 


where  K2  is  given  by  (4.6) 


(4.16a) 

(4.16b) 


K2  = -8*2X2  =  (4.16c) 

and  L2  is  given  by  (4.12) 

L2  =  -Y2c;  = 

Note  the  separation  structure  of  this  solution.  The  computation  of  K2  does  not  depend 
upon  the  influence  of  the  disturbances,  which  are  accounted  for  in  Bi  and  C2-  The 
computation  of  L2  does  not  depend  upon  the  weightings  in  the  cost  function,  which  are 
accounted  for  in  Ci,  or  the  manner  in  which  the  control  u  affects  the  state,  which  is 
accounted  for  in  ^2.  In  other  words,  the  problem  of  control  and  the  problem  of  state 
estimation  are  entirely  decoupled. 

Note  also  that,  in  order  to  arrive  at  the  (relatively)  simple  control  equations  described 
by  Doyle  et  al  (1989)  and  outlined  in  this  section  and  the  next,  the  matrices  A,  Bi,  B2-> 
Cl,  C2,  Di2j  and  D21  are  assumed  to  satisfy  eight  required  properties.  The  first  four  of 
these  properties 


{A^Bi)  stabilizable  (Ci,A)  detectable 

(^,^2)  stabilizable  {C2,A)  detectable 

are  verified  a  posterioriy  simply  by  examination  of  the  results.  (Note  that  the  analysis 
of  §3.3  indicates  that  there  is  only  one  slightly  unstable  mode  for  this  system,  and  that 
this  mode  is  both  sensitive  to  the  application  of  control  and  easily  discerned  by  the  mea¬ 
surements.  Thus,  we  may  presume,  but  not  assert  rigorously,  that  these  four  conditions 
will  in  fact  be  satisfied.)  The  matrices  are  constructed  to  satisfy  the  other  four  of  these 
properties  identically 


Cl  =  0  Di2  =  I 

Bid;,^o  ^21^2*1-/, 

as  may  be  verified  directly  with  the  definitions  of  these  matrices. 

4.2.  Hoo  control 

The  Hoo  estimator /controller  described  in  this  section  is  very  similar  to  the  H2  estima¬ 
tor/controller  described  previously.  Consideration  is  now  given  to  disturbances,  which 
we  shall  distinguish  with  a  new  variable  %,  of  the  “worst”  possible  structure  (as  made 
precise  below),  rather  than  the  Gaussian  structure  assumed  in  the  H2  case.  Considered 
in  the  frequency  domain,  the  estimator/controllers  developed  in  this  section  provide  a 
system  behaviour  in  which  the  maximum  singular  value  of  the  closed-loop  transfer  func¬ 
tion,  also  known  as  the  “oo-norm”,  is  less  than  some  constant,  which  shall  be  referred  to 
as  7.  As  this  approach  may  be  interpreted  as  bounding  the  oo-norm  of  the  transfer  func¬ 
tion  from  the  disturbances  to  the  performance  measure,  it  is  referred  to  as  Hoo  control. 
For  further  details  of  the  frequency- domain  explanation  of  Hooy  the  reader  is  referred  to 
Doyle  et  al.  (1989)  and  Zhou,  Doyle,  &  Glover  (1996). 
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The  governing  equations  to  be  considered  in  this  section  are  identical  to  (4.15): 

X  =  Ax  +  Bi  ^+^2  'fJ'  (4.17a) 

2  =  Ci.T+  Duu  (4.17b) 

y-(72.7;  +  T>2iX.  (4.17c) 


As  before,  the  Gi  and  G2  matrices  used  to  define  this  system  describe  any  covariance 
structure  of  the  disturbances  known  or  expected  a  priori  (for  instance,  if  one  measurement 
is  known  to  be  noisier  than  another).  These  matrices  are  taken  as  identity  matrices  if  no 
such  structure  is  known  in  advance. 

An  Hoo  estimator/controller  is  sought  to  relate  the  observations  y  to  the  control  li, 
which  is  applied  to  control  the  evolution  of  the  state  x  such  that  the  cost  Joo{z)  is 
minimized  in  the  presence  of  some  “worst  case”  disturbance  %.  The  flow  of  information 
is  represented  by  a  block  diagram  similar  to  that  shown  in  §4.1.3. 

Effectively,  the  cost  function  considered  for  Hoo  control  is 

Joo  =  E[X*  Qx  +  e  u*u  -  7"  x*x]-  (4.18) 

A  u  is  sought,  through  an  estimator/controller  of  the  form  (4.2) — (4.3),  to  rndnirndze 
Jooi  while  simultaneously  an  external  disturbance  x  is  sought  to  ma.ximize  (In  this 
manner,  %  is  the  “worst  possible”  disturbance,  as  it  is  exactly  that  disturbance  which 
increases  the  relevant  cost  function  the  most.)  Thus,  the  TLoo  problem  is  a  “min-max” 
problem.  The  term  involving  —7^  limits  the  magnitude  of  the  unstructured  disturbance 
in  the  maximization  of  Joo  with  respect  to  x  in  a  manner  analogous  to  the  term  involving 
,  which  limits  the  magnitude  of  the  control  in  the  minimization  of  Joo  with  respect  to 
u. 

The  result  (Doyle  et  ai  1989)  is  that  an  Hoo  estimator/controller  of  the  form  (4.2) — 
(4.3)  which  minimizes  Jo©  in  the  presence  of  some  component  of  the  worst  case  unstruc¬ 
tured  disturbance  x.  this  system  is  given  by 

X  =  AxA  B2U-  Loo  iv  ~  C2  x)  (4.19a) 

u  —  KoqX  (4.19b) 

where  Koo  is  given  by 

k^  =  -b;x^  = 

and  Loo  is  given  by 

Loo  =  -YooC;  - Ric  7"' Q Cl C2 j 

Note  first  that,  in  the  7  — >  cx)  limit,  the  H2  estimator/cont roller  is  recovered,  so  the  set  of 
two  Riccati  equations  in  (4.19)  describes  both  the  H2  (optimal  control  -f  Kalman-Bucy 
filter)  and  the  Hoo  problems. 

It  may  also  be  shown  that,  as  the  upper-right  blocks  of  the  Hamiltonians  may  not  be 
negative  definite,  a  solution  to  these  Riccati  problems  exists  only  for  sufficiently  large  7; 
the  smallest  7  =  70  for  which  a  solution  to  these  equations  exists  may  be  found  by  trial 
and  error  (Doyle  et  al  1989).  An  Hoo  estimator /controller  for  7  >  70  is  referred  to  as 
suboptimal. 
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4.3.  Com>parison  of  H2  and  Hoc  control  equations 
Most  of  the  robustness  problems  associated  with  H2  stem  from  the  state  estimation.  Op¬ 
timal  (LQR)  controllers  themselves,  provided  with  full  state  information,  generally  have 
excellent  performance  and  robustness  properties  (Dailey  et  al.  1990).  Recall  from  §4.1.3 
that  the  problems  of  control  and  state  estimation  in  the  %  formulation  are  decoupled. 

An  important  observation  of  §4.2  is  that  the  problems  of  control  and  state  estimation 
in  the  formulation  are  coupled.  Specifically,  the  computation  of  depends  on  the 
expected  covariance  of  the  state  disturbances,  which  are  accounted  for  in  Ri,  and  the 
computation  of  Loo  depends  on  the  weightings  in  the  cost  function,  which  are  accounted 
for  in  Cl ,  This  is  one  of  the  essential  features  of  ?{oo  control. 

By  taking  into  account  the  expected  covariance  of  the  state  disturbances,  reflected  in 
Ri,  when  determining  the  state  feedback  matrix  Koa  the  components  of  x  correspond¬ 
ing  to  the  components  of  x  that  are  expected  to  have  the  smallest  forcing  by  external 
disturbances  are  weighted  least  in  the  feedback  control  relationship  u  =  Koo^- 

Similarly,  by  taking  into  account  the  weightings  in  the  cost  function,  reflected  in  Ci , 
when  determining  the  estimator  feedback  matrix  Loo  ?  the  components  of  x  corresponding 
to  the  components  of  x  that  are  least  important  in  the  computation  of  Joo  are  forced 
with  the  smallest  corrections  by  the  output  injection  term  Loo  {y  —  y)  in  the  equation  for 
the  estimator. 

By  applying  strong  control  only  on  those  components  of  x  significantly  excited  by 
external  disturbances,  and  by  applying  strong  estimator  corrections  only  to  those  com¬ 
ponents  of  X  important  in  the  computation  of  the  cost  function,  7ioo  feedback  gains  for 
components  of  the  system  not  relevant  to  the  control  problem  are  reduced  from  those 
in  the  ^2  case.  With  such  feedback  gains  reduced,  the  stability  properties  of  Hoo  es¬ 
timator/controllers  in  the  presence  of  state  disturbances  and  measurement  noise  may 
be  expected  to  be  better  than  their  %2  counterparts,  at  the  cost  of  a  (hopefully,  small) 
degradation  of  performance  in  terms  of  the  2-norm  of  the  output  ^  for  the  undisturbed 
system. 

4.4.  Numerical  method 

Standard  numerical  techniques  are  now  applied  to  all  aspects  of  this  problem.  In  order 
to  simplify  both  the  theory  to  be  presented  and  the  numerical  algorithm  to  be  coded,  no 
further  manipulation  of  the  equations  is  used  beyond  the  matrix  representations  (4.17) 
and  (4.19).  It  was  observed  that  the  minimal  realization  approach  (Kailath  1980)  is  well 
suited  to  reduce  the  computation  time  necessary  to  determine  eflFective  control  algorithms 
by  the  present  approach;  however,  such  an  approach  was  not  found  to  be  necessary  in 
the  present  case. 

The  algebraic  Riccati  equations  are  solved  using  the  method  of  Laub  (1991),  which 
involves  a  Schur  factorization.  This  is  found  to  be  a  stable  numerical  algorithm  for  all 
cases  tested.  The  implementation  of  Laub’s  method  is  written  in  Fortran-90  and  follows 
closely  the  algorithm  used  by  the  Matlab  function  are  .m  (Grace  et  al.  1992).  A  Lyapunov 
solver,  modelled  after  the  Matlab  function  lyap.m,  is  also  used  to  compute  the  system 
Gramians. 

Two  LAPACK  routines  (Anderson  et  al.  1995),  zgeev.f  and  zgees.f,  are  used  to 
compute  eigenvalues/eigenvectors  and  Schur  factorizations.  These  routines  are  compiled 
in  quad  precision  (128  bits  per  real  number)  to  ensure  sufficient  numerical  precision  in  the 
eigenvalue  computation.  All  computations  are  carried  out  with  W  =  140  to  ensure  good 
resolution  of  all  significant  eigenmodes.  The  eigenvalues  of  A  match  all  those  tabulated 
by  Orszag  (1971)  to  all  eight  decimal  places,  as  shown  in  table  1,  indicating  that  this 
numerical  method  is  sufficiently  accurate. 
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5.  Performance  of  controlled  systems 

We  now  examine  the  behaviour  of  the  “closed-loop”  systems  obtained  by  application 
of  the  above  controllers  and  estimators  to  the  “nominal”  {i.e.  no  disturbances)  channel 
flow  stability  problem.  In  other  words,  we  examine  the  behaviour  of  the  flow  and  the 
estimator/cont rollers  operating  together  as  a  single  dynamical  system.  By  looking  at 
“root  locus”  plots  which  map  the  movement  of  the  eigenvalues  of  these  systems  in  the 
complex  plane  with  respect  to  the  relevant  parameters,  this  behaviour  is  well  quantified, 
as  long  as  all  perturbations  remain  small  enough  that  the  linearity  assumption  remains 
valid.  We  shall  also  examine  the  control  and  observation  sensitivities  defined  in  §3.2  for 
two  special  cases  in  order  to  better  understand  the  fundamental  limitations  of  controllers 
and  estimators  applied  to  the  present  system. 


5.1.  712  control 

5.1.1.  Optim.al  control  (LQR) 

In  order  to  investigate  the  controllability  of  the  closed- loop  eigenmodes  when  all  modes 
are  observable,  consider  the  system  described  in  §4.1.1.  With  r  =  0  and  examining  only 
the  equations  for  x  and  .t,  the  plant  is  given  (in  the  shorthand  notation  used  in  §4)  by 


Vlqr  = 


X 

u 

X 

■  A 

B2 

X 

I 

0 

with  the  control  now  given  by 


u  —  K2  X  -|-  7/, 


where  an  additional  control  term  u'  has  been  added  to  study  the  sensitivity  of  the  closed- 
loop  system  to  further  modification  of  the  control.  Putting  the  plant  and  the  controller 
together,  the  closed-loop  system  may  be  represented  by 


LQR  (closed  loop) 


X 

u' 

.7; 

A  +  B2  K2 

B2 

X 

I 

0 

(5.1) 


The  eigenmodes  of  A/Ca  A  +  B2  K2  describe  the  dynamics  of  the  closed-loop  system  for 
the  unmodified  control  rule  (w'  =  0).  Figure  2a  shows  the  movement  of  these  eigenvalues 
with  respect  to  the  free  parameter  of  the  control  problem,  £,  used  to  determine  K2>  The 
eigenvalues  for  ^  00  are  observed  to  be  very  near  those  of  the  uncontrolled  system  A 

in  figure  la,  with  the  previously  unstable  mode  moved  just  to  the  left  of  the  imaginary 
axis.  The  eigenvalues  generally  move  to  the  left  as  i  is  decreased.  Figure  2b  shows  the 
shape  of  the  first  four  eigenmodes  of  the  closed-loop  system.  Comparing  figure  2b  with 
figure  lb,  it  is  seen  that  the  control  modifies  most  those  eigenmodes  with  significant 
variations  near  the  wall. 

The  sensitivity  of  the  eigenmodes  of  the  system  (5.1)  to  modification  of  the  control 
rule  may  be  quantified  by  performing  the  analysis  of  §3.2.1,  replacing  the  eigenmodes  of 
A  by  the  eigenmodes  of  Ak2-  The  result  of  this  analysis  for  small  £  is  shown  in  table  2. 
This  table  shows  that,  in  the  £  — )•  0  limit,  the  system  matrix  is  modified  to  the  point  that 
the  eigenmodes  are  no  longer  sensitive  to  further  modification  of  the  control.  In  other 
words,  all  the  controllable  dynamics  of  the  system  have  been  modified  by  K2  and  are 
accounted  for  in  the  closed  loop  system  in  this  limit.  This  is  one  demonstration  that  the 
optimal  controller  extracts  the  best  possible  performance  from  a  given  (full-information) 
system. 
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3 

A, 

/i 

3 

-0.03513233  - 

0.96462128z' 

0.000000029 

4 

-0.03518652  - 

0.96464261Z 

0.000000001 

5 

-0.06255259  - 

0.29262711* 

0.000001101 

6 

-0.06310358  - 

0.93629329i 

0.000000070 

7 

-0.06325089  - 

0.93635257i 

0.000000003 

1 

-0.06644730  - 

0.29721403* 

0.000001116 

8 

-0.09102975  - 

0.9079395H 

0.000000129 

9 

-0.09130964  - 

0.90805917* 

0.000000008 

10 

-0.11890731  - 

0.87955083* 

0.000000226 

11 

-0.11936036  - 

0.87976246* 

0.000000020 

12 

-0.14335180  - 

0.43962023* 

0.000002303 

14 

-0.14673294  - 

0.85111508* 

0.000000414 

15 

-0.14739907  - 

0.8514616H 

0.000000045 

13 

-0.14803996  - 

0.44586838* 

0.000003081 

16 

-0.17450455  - 

0.82261690*' 

0.000000842 

Table  2.  Least  stable  eigenmodes  of  the  closed-loop  system  Ak2  and  their  sensitivity  to  control 
for  the  optimal  controller  in  the  cheap  control  limit  {£  =  10“^).  The  numbering  of  the  eigenvalues 
shown  is  the  same  as  the  numbering  of  the  eigenvalues  of  table  1  to  which  they  are  connected  by 
the  root  locus  of  figure  2.  Note  that  the  control  in  this  limit  drives  all  eigenmodes  to  positions 
at  which  they  are  insensitive  to  further  modifications  of  the  control,  as  illustrated  by  the  large 
reductions  in  fj.  Note  also  that  those  eigenmodes  with  the  largest  values  of  fj  in  table  1 
(specifically,  those  in  the  lower  branch)  have  moved  the  most. 


5.1.2.  Kalman-Bucy  filter  (KBF) 

The  estimator  itself  has  its  own  set  of  dynamics.  These  dynamics  are  captured  by 
the  equations  for  the  state  estimator  error,  as  described  in  §4.1.2.  We  now  make  use  of 
this  system  in  order  to  investigate  the  observability  of  closed-loop  eigenmodes  when  all 
modes  are  controllable.  With  tt;  =  0  and  examining  only  the  equations  for  xe  and 
this  plant  is  given  by 


“Pkef  = 

with  the  output  injection  now  given  by 


XE 

u 

XE 

■  A 

I 

VE 

C2 

0 

u  =  L2  VE  +  u  , 


where  an  additional  output  injection  term  u'  has  been  added  to  study  the  sensitivity  of 
the  closed-loop  system  to  further  modification  of  the  output  injection  rule.  Putting  the 
plant  and  the  estimator  together,  the  closed-loop  system  may  be  represented  by 


KBF  (closed  loop)  — 


XE 

u 

XE 

A  -{-  Z/2  C2 

I 

VE 

C2 

0 

(5.2) 


The  eigenmodes  of  ^^2  =  A L2  C2  describe  the  dynamics  of  the  closed-loop  system  for 
the  unmodified  output  injection  rule  {u'  —  0).  Figure  3  shows  the  movement  of  these 
eigenvalues  with  respect  to  the  free  parameters  of  the  estimator  problem.  (This  is  done 
by  assuming  that  the  matrices  describing  the  covariance  of  the  disturbances  have  the 
simple  form  Gi  =  gi  I  and  G2  —  g2  I ^  where  <71  and  g2  are  real  scalars.)  The  eigenvalues 
for  =  ^2  — ^  0  are  very  near  those  of  the  uncontrolled  system  A  in  figure  la,  with  the 
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Figure  3.  Root  locus  of  least  stable  eigenvalues  of  Al^  as  a  function  of  the  free  parameters 
of  the  7/2  estimator,  (ji  and  <72  (note  that  we  take  gi  g2  for  the  purpose  of  drawing  the  root 
locus).  The  eigenvalues  for  =  ^^2  ^  0  are  marked  with  an  (x). 
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8 
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0.000000673 

9 
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0.000000011 

1 
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0.000000094 

10 
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0.000000646 

11 

-0.11936807  - 

0.879757092 

0.000000014 

12 

-0.14209547  - 

0.25910275t 

0.000000130 

14 

-0.14584717- 
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15 

-0.14741926  - 
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0.000000014 

16 

-0.17347707  - 
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0.000000399 

13 

-0.17418920  - 

0.403146561 

0.000002002 

Table  3.  Least  stable  eigenmodes  of  the  closed-loop  system  Al2  their  sensitivity  to  obser¬ 
vation  for  the  Kalman-Bucy  filter  in  the  large  disturbance  limit  (^1  =  ^2  =  10^).  The  numbering 
of  the  eigenvalues  shown  is  the  same  as  the  numbering  of  the  eigenvalues  of  table  1  to  which 
they  are  connected  by  the  root  locus  of  figure  1.  Note  that  the  estimator  in  this  limit  modifies 
all  eigenmodes  until  the  measurements  are  no  longer  sensitive  to  them,  as  illustrated  by  the 
large  reductions  in  gj.  Note  also  that  those  eigenmodes  with  the  largest  values  of  gj  in  table  1 
(specifically,  those  in  the  lower  branch)  have  moved  the  most. 
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Figure  4.  Least  stable  eigenvalues  of  the  composite  closed-loop  system  with  the  %  estima¬ 
tor/controller,  taking  I  =  gi  z=  g2  =  1,  Note  that  the  eigenvalues  are  simply  the  eigenvalues  of 
the  closed  loop  controller  (-{-)  together  with  those  of  the  closed  loop  estimator  (*). 

previously  unstable  mode  moved  just  to  the  left  of  the  imaginary  axis.  The  eigenvalues 
generally  move  to  the  left  as  gi  and  g2  are  increased. 

The  sensitivity  of  measurements  ije  the  the  eigenmodes  of  the  system  (5.2)  may 
be  quantified  by  performing  the  analysis  of  §3.2.2,  replacing  the  eigenmodes  of  A  by  the 
eigenmodes  of  Al^.  The  result  of  this  analysis  for  large  =  g2  is  shown  in  table  3. 
This  table  shows  that,  in  the  gi  =  g2  ^  oo  limit,  the  system  matrix  is  modified  to  the 
point  that  the  measurements  are  no  longer  sensitive  to  the  eigenmodes  of  the  closed- 
loop  system.  In  other  words,  all  the  measurable  dynamics  of  the  system  have  been 
extracted  by  L2  and  are  accounted  for  in  the  closed  loop  system  in  this  limit.  This  is 
one  demonstration  that  the  Kalman-Bucy  filter  extracts  the  best  possible  state  estimate 
from  a  given  (fully-controllable)  state  estimator. 

5.1.3.  7^2  control  (LQG  =  LQR -h  KBF) 

It  was  mentioned  in  §4.1.3  that  the  estimator/ controller  which  minimized  the  relevant 
cost  functional  {J2)  in  the  presence  of  Gaussian  disturbances  could  be  found  by  consid¬ 
ering  the  controller  and  estimator  problems  separately.  In  this  section,  it  is  shown  that 
the  closed-loop  performance  of  a  system  of  the  form  (4.15)  (without  disturbances) 

X  =  A  .7;  +  Z?2 
y  =  C2X 

combined  with  an  estimator/controller  of  the  form  (4.16) 

X  =  Ax B2  u  ~  L2  [y  -  C2  x) 
u  =  K2  X 

may  also  be  evaluated  by  considering  the  estimator  and  controller  problems  separately. 
To  accomplish  this,  simply  combine  the  above  equations  into  the  closed-loop  composite 
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Figure  5.  Root  locus  of  least  stable  eigenvalues  of  the  Roo  controller  versus  7,  taking  I  =  100, 
g-^  =  g2  =  0.001.  The  result  with  7  ^  00,  marked  with  the  (x),  gives  the  corresponding 
controller.  Note  that  the  Roo  controller  modifies  only  the  least  stable  eigenmode  of  this  R2 
result,  without  expending  any  extra  control  effort  to  control  those  eigenmodes  not  associated 
with  the  maximally  unstable  component  of  the  system.  Note  also  that  7  =  70,  marked  with 
the  (o),  is  reached  by  reducing  7  until  the  least  stable  eigenvalue  corresponds  to  one  of  the 
uncontrollable  eigenmodes  in  the  upper  branch,  which  cannot  be  moved  further  left;  in  the 
present  case,  this  corresponds  to  a  numerical  value  of  70  =  0.26. 

system 

A  B2  K2 

—L2  O2  ^  B2  K2  +  L2  O2 

Gaussian  elimination,  first  on  the  rows  and  then  on  the  columns,  reveals  that  the  eigen¬ 
values  of  this  system  are  the  same  as  the  eigenvalues  of  the  system 

^  +  B2  K2  B2  K2  \ 

0  A  L2  C2J 

In  other  words,  the  eigenvalues  of  the  closed-loop  composite  system  for  the  R2  problem 
are  simply  the  union  of  the  eigenvalues  of  the  controlled  system  =  A B2  K2  and 
the  eigenvalues  of  the  estimated  system  Al^  A  A  L2C2  discussed  in  the  previous  two 
sections  and  illustrated  in  figure  4. 

5.2.  Roo  control 

As  with  the  R2  estimator/controller,  the  performance  of  the  closed  loop  composite  system 
with  the  Roo  estimator/controller 

A  B2  Koc 

—Loo  C2  A-\-  B2  Koo  +  Loo  C2 

may  be  evaluated  by  considering  the  performance  of  the  controlled  system  Ajc^  =  A -\- 
B2  Koo  and  the  performance  of  the  estimated  system  Ai^  =  A  -(-  Loo  C2  separately. 
The  root  locus  of  the  eigenvalues  of  Ak^  are  plotted  with  respect  to  the  parameter  7 
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of  the  Hoo  problem  in  figure  5,  clearly  illustrating  the  tendency  of  %oo  controllers  to 
modify  only  the  least  stable  components  of  the  system,  as  opposed  to  the  %2  controller 
of  figure  2,  which  modifies  all  controllable  modes  of  the  system. 


6.  Conclusions 

Optimal  and  robust  control  theories  have  been  successfully  applied  to  the  Orr-Sommerfeld 
equation.  Given  control  on  the  wall-normal  component  of  boundary  velocity  only,  the 
flow  system  is  shown  to  be  stabilizable  but  not  controllable.  Given  measurements  of  wall 
skin-friction  only,  the  flow  system  is  shown  to  be  detectable  but  not  observable.  It  is 
shown  that  TL2  controllers /estimators  modify  all  of  the  controllable/observable  modes 
of  the  system.  In  contrast,  the  7/oo  controllers  modify  the  corresponding  ?^2  controllers 
only  in  the  most  unstable  component,  as  Hoc  targets  a  bound  only  on  the  maximum 
value  of  the  transfer  function. 

In  the  £  0  limit  of  the  %2  controller,  corresponding  to  cheap  control  and  thus  large 

values  of  w,  all  eigenmodes  of  the  closed-loop  controlled  system  are  shown  to  be  modified 
to  points  at  which  they  are  no  longer  sensitive  to  further  modifications  of  the  control. 
Similarly,  in  the  g\  g2  00  limit  of  the  7^2  estimator,  accounting  for  large  disturbances 
on  both  the  state  and  the  measurements,  all  eigenmodes  of  the  closed-loop  system  for 
the  estimator  error  are  shown  to  be  modified  to  points  at  which  they  are  not  discernible 
by  flow  measurements. 

These  results  indicate  that  'H2  controllers  and  estimators  are  optimal  for  their  desired 
purposes,  but  may  contain  large  feedback  gains.  On  the  other  hand,  7£oo  controllers 
only  target  the  least  stable  components  of  the  system,  and  thus  have  smaller  feedback 
gains  while  still  achieving  the  same  worst  case  performance  for  the  nominal  plant.  Such 
reduced  feedback  gains  generally  result  in  improved  robustness  to  inaccuracies  in  the 
system  model. 
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PART  B. 

Optimal  control  of  turbulence 

Optimal  control  theory  is  used  to  determine  controls  that  effectively  reduce  the  drag  of 
a  turbulent  flow  in  a  plane  channel.  Wall  transpiration  (unsteady  blowing/suction)  with 
zero  net  mass  flux  is  used  as  the  control.  The  technique  described  is  unique  from  the 
standpoint  that  it  is  mathematically  based  solely  on  the  control  objective,  the  equations 
governing  the  fluid  flow,  and  instantaneous  observations  of  the  flow,  without  the  ad  hoc 
procedures  normally  used  to  accomplish  control  of  complex  nonlinear  phenomena  such 
as  turbulence.  Drag  reduction  of  over  50  percent  is  obtained  using  an  optimal  controller 
in  a  direct  numerical  simulation  of  a  turbulent  channel  flow  at  Rcr  =  180,  which  far 
exceeds  what  has  been  obtained  to  date  via  adaptive  and  intuition-based  control  rules  in 
similar  flows. 

The  algorithm  used  is  computationally  intensive  and  requires  full  flowfleld  information, 
and  therefore  can  not  be  implemented  in  a  laboratory.  However,  these  calculations  allow 
us  to  quantify  the  best  possible  system  performance  given  a  certain  class  of  flow  actuation 
and  qualitatively  identify  how  optimized  controllers  interact  with  the  coherent  structures 
of  the  turbulence.  In  so  doing,  an  important  step  is  made  in  the  progression  towards 
practical  and  efficient  turbulence  control  strategies  based  on  optimal  control  techniques. 


1.  Introduction 

The  recent  development  of  the  technology  necessary  to  produce  micro-scale  mechan- 
canical  devicesf,  commonly  referred  to  as  Microfabricated  Electro-Mechanical  Systems 
(MEMS),  has  prompted  researchers  to  revisit  questions  heretofore  thought  to  be  purely 
academic.  The  present  question  is  exactly  of  this  nature:  assuming  that  individual  small- 
scale  turbulent  fluctuations  may  somehow  be  measured  and  that  concomitant  small-scale 
forcing  of  the  turbulent  fluid  may  somehow  be  attained,  how  much  do  practical  engineer¬ 
ing  designs  stand  to  benefit,  where  should  the  control  be  applied,  and  what  control 
algorithms  are  most  effective?  The  present  work  attempts  to  cast  these  questions  in  a 
rigorous  framework,  present  a  mathematical  approach  for  their  solution,  and  demonstrate 
an  effective  implementation  of  the  control  algorithm  in  a  fully-developed  turbulent  flow. 

We  first  briefly  summarize  recent  approaches  to  determine  effective  turbulence  control 
algorithms,  categorizing  these  approaches  to  the  feedback  control  problem  by  examining 
their  mathematical  dependence  on  the  equations  governing  the  system.  This  discussion 
puts  the  present  approach  in  context  with  other  techniques  currently  under  investigation. 
For  a  more  thorough  discussion  along  this  line,  see  Moin  &  Bewley  (1994). 

1.1.  Adaptive  netv)orks 

The  first  class  of  schemes  which  may  be  proposed  to  achieve  small  scale  flow  control 
actually  makes  no  explicit  reference  to  the  dynamics  known  to  take  place  in  the  flow  or 
the  Navier-Stokes  equations  known  to  govern  these  dynamics.  Instead,  a  “reasonable” 
network  is  fashioned  which  takes  as  input  those  measureable  flow  quantities  assumed  to 

f  For  a  recent  review  of  this  subject  as  it  applies  to  fluid  mechanics,  see  McMichael  (1996). 
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be  most  relevant  to  the  control  problem  and  produces  as  output  the  requisite  control 
velocity.  The  coefficients  of  this  network  are  then  “trained”  by  applying  the  control 
network  to  the  flow  and  gradually  adjusting  the  coefficients  in  a  heuristic  manner  based 
on  the  resulting  evolution  of  the  flow.| 

As  an  example  of  one  adaptive  approach,  an  adaptive  inverse  technique  has  been 
applied  by  Lee  et  al  (1996)  to  a  turbulent  channel  flow  at  Rcr  =  100,  providing  ap¬ 
proximately  20  percent  drag  reduction.  This  approach  first  develops  an  approximate 
“inverse”  model  between  measurable  flow  quantities  (as  input)  and  the  control  forcing 
(as  output)  with  an  adaptive  technique.  Each  iteration  of  the  adaptation  for  this  inverse 
model  consists  of  three  steps:  1)  computing  the  error  of  the  model  output  with  respect  to 
the  desired  model  output  (the  actual  control  forcing  used),  2)  determining  the  influence 
of  the  various  weights  in  the  model  on  this  error,  then  3)  updating  all  the  weights  in  the 
model  a  small  amount  in  a  manner  that  reduces  the  error.  When  applied  to  the  nonlinear 
adaptive  networks  commonly  used  for  this  purpose,  known  as  “neural  networks” ,  this  is 
commonly  referred  to  as  “back- propagation”  of  the  error.  Once  the  approximate  inverse 
model  between  the  flow  measurements  and  the  control  converges,  the  inverse  model  is 
used  to  compute  a  control  which  will  drive  the  flow  measurements  to  some  desired  state. 
In  the  case  of  Lee  et  al.  (1996),  the  desired  state  is  chosen  to  be  a  state  with  reduced 
spanwise  drag  fluctuations,  and  the  inverse  model  is  continually  trained  as  the  flow  sys¬ 
tem  evolves  in  time.  Such  on-line  training  of  the  controller  helps  to  provide  “robustness” 
to  possible  changes  in  the  dynamics  of  the  flow  system  to  be  controlled. 


1.2.  Intuition-based  approaches 

In  situations  in  which  the  dominant  physics  is  well  understood,  judgment  can  guide 
an  engineer  to  design  effective  control  schemes.  Success  is  limited,  however,  by  the 
engineer’s  understanding  of  the  physical  processes  involved;  in  the  case  of  turbulence, 
our  understanding  is  still  limited  despite  several  decades  of  intense  research. 

An  active  cancellation  scheme  was  used  by  Choi,  Moin,  &  Kim  (1994),  to  reduce  the 
drag  in  a  fully- developed  turbulent  flow  by  mitigating  the  effect  of  the  near- wall  vortices. 
By  opposing  the  near- wall  motions  of  the  fluid,  which  are  caused  by  the  near- wall  vortices, 
with  an  opposing  wall  control,  the  high  shear  region  was  lifted  away  from  the  wall.  A 
direct  numerical  simulation  of  this  scheme  applied  to  turbulent  channel  flow  at  Re^  =  100 
demonstrated  23  percent  drag  reduction  when  the  control  was  chosen  to  oppose  the 
vertical  motion  at  y'^  =  10. 

Sensing  the  instantaneous  normal  velocity  at  y~^  =  10  is,  of  course,  very  difficult. 
From  a  practical  standpoint,  it  is  highly  desirable  to  confine  both  sensing  and  actuation 
to  the  wall.  Thus,  Choi,  Moin,  &  Kim  (1994)  computed  the  correlation  of  quantities 
measurable  at  the  wall  with  the  normal  velocity  above  the  wall.  Surprisingly,  the  wall 
pressure  did  not  exhibit  a  high  correlation  with  the  normal  velocity.  Using  a  Taylor 
series  expansion  and  the  equation  of  continuity,  they  obtained  an  expression  relating  the 
normal  velocity  at  a  point  near  the  wall  to  the  instantaneous  wall  shear.  However,  using 
this  expression  to  estimate  the  normal  velocity  away  from  the  wall  resulted  in  only  a  6 
percent  drag  reduction.  This  is  comparable  to  the  drag  reduction  that  can  be  achieved 
with  simpler,  passive  means  such  as  riblets. 

f  Note  that  the  network  used  for  this  purpose  can  take  any  of  several  linear  or  nonlinear  forms. 
Hertz,  Krogh,  &  Palmer  (1991)  contains  a  survey  of  these  networks  and  outlines  strategies  for 
their  training. 
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1.3.  Dynamical  systems 

The  tools  of  dynamical  systems  theory  have  proven  useful  in  analyzing  and  interpreting 
turbulence  dynamics  (Aubry  et  al.  1988).  Due  to  their  large  range  of  spatial  and  tempo¬ 
ral  scales,  turbulent  flows  are  known  to  have  relatively  high  dimensions  in  this  framework 
even  at  fairly  low  Reynolds  numbers,  which  makes  analysis  of  these  systems  quite  difficult 
(Keefe,  Moin,  &  Kim  1992).  However,  there  has  been  some  instructive  work  in  represent¬ 
ing  the  dynamics  of  coherent  structures  in  wall-bounded  flows  with  systems  of  much  lower 
dimension  using  the  proper  orthogonal  decomposition  (POD)  method  (Berkooz,  Holmes, 
&  Lumley  1993).  This  method  provides  a  (numerically  determined)  set  of  eigenfunctions 
which  are  particularly  efficient  in  representing  second  order  turbulence  statistics  with  a 
small  number  of  modes  (perhaps  as  few  as  10  to  20)  in  the  cross-flow  plane. 

In  the  dynamical  systems  framework,  the  movement  of  the  near-wall  longitudinal  vor¬ 
tices  when  observed  in  a  cross-flow  plane  may  be  represented  locally  as  the  orbiting  of 
a  low  dimensional  state  around  several  unstable  fixed  points;  the  passage  of  one  set  of 
coherent  structures  leads  to  a  rapid  jump  in  the  state  to  a  different  unstable  orbit,  or  to 
a  different  distribution  of  near- wall  longitudinal  vortices  in  the  cross-flow  plane.  Such  a 
rapid  transition  between  critical  points,  followed  by  a  quiescent  period  in  which  the  flow 
pattern  remains  largely  unchanged,  is  referred  to  as  a  heteroclinic  cycle. 

Coller,  Holmes,  &  Lumley  (1994a,b)  consider  the  control  of  an  interesting  model  prob¬ 
lem  governed  by  a  simple  two-component  equation  with  similar  dynamics  to  this  model 
of  near- wall  longitudinal  vortices  {i.e.  attracting  heteroclinic  cycles)  subject  to  random 
excitation  to  account  for  unmodelled  disturbances.  They  develop  and  demonstrate  a 
strategy  which  delays  heteroclinic  transitions  in  this  model  as  long  as  possible  by  sensing 
when  the  state  is  near  an  unstable  fixed  point  and  maintaining  it  there  with  feedback 
control  for  as  long  as  possible.  Once  the  state  diverges  from  this  fixed  point,  presum¬ 
ably  due  to  the  unmodelled  dynamics  of  the  flow,  control  is  turned  off  until  the  state 
approaches  the  neighborhood  of  another  unstable  fixed  point. 


1.4.  Rigorous  optimization  of  practical  control  algorithms — a  preview 

The  present  work  is  a  first  step  towards  developing  practical  controllers  of  the  two  types 
illustrated  in  figure  1.  The  (initially  undetermined)  coefficients  present  in  both  con¬ 
figurations  may  be  optimized  rigorously  with  approaches  based  on  the  adjoint  analysis 
developed  in  this  paper,  and  are  still  under  development  (Bewley,  Moin,  &  Temam  1996). 

In  the  output  feedback  configuration,  the  flow  is  controlled  using  computationally 
inexpensive  direct  feedback  from  instantaneous  flow  measurements.  The  structure  of  the 
controller  may  be  nonlinear  and  may  incorporate  a  finite  impluse  response  (FIR)  filter 
to  account  for  information  from  past  measurements  in  the  control  rule. 

In  the  estimator /controller  configuration,  a  time-evolving  estimate  of  the  flow  state 
near  the  wall  is  first  developed,  effectively  accumulating  the  information  reflected  by  the 
stream  of  measurements  from  a  few  noisy  sensors.  The  flow  is  then  controlled  with  a 
(possibly  nonlinear)  control  rule  based  on  this  flow  estimate.  The  estimation  problem 
and  the  control  problem  become  linked  when  they  are  optimized  in  the  presence  of  a 
small  component  of  “worst  case”  noise  which  maximally  aggrevates  the  coupled  system. 
Such  an  approach  is  well  developed  for  linear  problems,  and  is  referred  to  as  control; 
Doyle  et  al.  (1989)  presents  a  compact  form  of  this  approach  which  makes  it  straight¬ 
forward  to  apply  to  linear  problems,  as  illustrated  in  Bewley,  Agarwal,  &  Liu  (1996)  for 
the  control  of  the  linear  stages  of  transition.  Methods  to  extend  this  “robust”  approach 
to  nonlinear  problems  are  still  under  analysis. 
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a)  Output  feedback  configuration  b)  Estimator/ controller  configuaration 

Figure  1.  Examples  of  practical  control  configurations  based  on  a  few  noisy  flow  measurements. 
Both  closed-loop  configurations  may  be  optimized  with  techniques  based  on  the  approach  de¬ 
veloped  in  this  paper,  and  are  still  under  development. 


2.  Governing  equations 

The  flow  problem  we  will  consider  is  fully- developed  turbulent  channel  flow  with  no-slip 
walls  and  wall-normal  velocity  boundary  conditions  ^  applied  as  the  control.  Though 
this  is  an  idealized  geometry,  it  will  give  insight  into  the  turbulent  behavior  which  can 
later  be  exploited  in  more  practical  configurations,  such  as  the  control  of  a  spatially 
developing  boundary  layer  with  discrete  wall-mounted  actuators.  The  present  problem 
is  governed  by  the  unsteady,  incompressible  Navier-Stokes  equation  and  the  continuity 
equation  inside  the  domain  Q,  and  wall-normal  velocity  boundary  conditions  on  the  walls 
On  the  remainder  of  the  boundary  of  the  three-dimensional  volume  denoted 
5^2,  periodic  boundary  conditions  are  applied.  The  extent  of  the  computational  domain 
is  chosen  to  be  large  enough  in  the  wall-parallel  directions  that  the  convenient  (though 
non-physical)  periodic  boundary  conditions  on  80.2  have  no  effect  on  the  nature  of  the 
near- wall  turbulence,  as  illustrated  qualitatively  in  figure  2. 

The  governing  equations  may  be  written  functionally  as 

ff{U)  =  0  inO  (2.1a) 

with  boundary  conditions 

Ui  =  ^rii  on  (2.1b) 


and  prescribed  initial  conditions 


Ui  =  Ui{0)  at  ^  =  0. 


(2.1c) 


For  clarity,  all  differential  equations  will  be  written  in  operator  form  in  this  discussion, 
with  these  operators  defined  when  first  introduced.  The  (nonlinear)  Navier-Stokes  op¬ 
erator  for  the  present  case,  in  which  the  flow  is  assumed  to  have  uniform  density  and 
viscosity,  is  given  by 


mu)  = 


( dui  dui  d^Ui  dp  I  \ 

at  dxj  axj  p  oxi  p 

duj 

dxj 


35 


PART  B.  Optimal  control  of  Uirbulence 

In  this  discussion,  xi  is  the  streamwise  direction,  X2  is  the  wall-normal  direction,  x^ 
is  the  spanwise  direction,  the  Ui^s  are  the  corresponding  velocities,  p  is  the  pressure, 
p  is  the  density,  p  is  the  absolute  viscocity,  and  i/  =  p/p  is  the  kinematic  viscocity. 
The  flow  is  sustained  by  a  mean  pressure  gradient  Px  in  the  streamise  direction,  which  is 
modified  at  each  time  step  in  order  to  maintain  a  constant  mass  flux  through  the  channel. 
Define  also  n  as  a  wall-normal  unit  vector  directed  into  the  channel,  Tyj  =  P' dui / dn\^all 
as  the  mean  skin  friction  on  the  wall  for  the  uncontrolled  channel  (averaged  in  space 
and  time),  Ur  =  as  the  mean  friction  velocity,  S  as  the  channel  half-width, 

and  Rcr  =  UrSlj/  as  the  Reynolds  number  based  on  the  mean  friction  velocity  and  the 
channel  half  width.  The  flow  considered  in  this  work  is  taken  at  Rcr  =  180.  All  velocities 
are  normalized  by  the  friction  velocity  and  therefore  may  also  be  marked  with  a  (‘^) 
superscript.  All  lengths  are  normalized  by  5  unless  marked  with  a  (■^)  superscript,  in 
which  case  they  are  normalized  by  the  wall  unit  vjur^  All  times  are  normalized  by 
unless  marked  with  a  ("*')  superscript,  in  which  case  they  are  normalized  by  Note 

that,  with  this  normalization,  v  =  l/Rcr  in  the  above  equation  for  Af{U), 

Three  state  vectors  are  used  in  this  work:  the  flow  U,  the  flow  perturbation  U\  and 
the  adjoint  ?7*: 

^  _  /u'i{xi,X2,X3,t)\  _  /u’f{xi,X2,X3,t) 

\pixi,X2,X3,t)  J  '  \P' ixi,X2,X3,t)J  '  \p*  {xi,X2,X3,t) 

Note  that  each  of  these  vector  fields  is  composed  of  three  velocity  components  and  a  pres¬ 
sure  component.  The  motivation  for  introducing  C/'  and  U*  will  be  apparent  in  the  deriva¬ 
tion  of  the  control  equations  to  follow,  and  the  parital  differential  equations  governing 
these  fields  will  be  derived.  Only  after  the  control  problem  has  been  derived  in  diflFerential 
form  is  it  discretized  in  space  and  time.  For  the  current  three-dimensional  nonlinear  prob¬ 
lem,  this  approach  is  found  to  yield  systems  of  equations  which  are  easiest  to  understand 
and  to  code.  Note  that  for  simpler  systems  of  equations,  such  as  the  one-dimensional 
linear  problem  of  transition  control  examined  in  Bewley,  Agarwal,  &  Liu  (1996),  the  ma¬ 
trix  control  equations  derived  from  the  discrete  form  of  the  governing  equations  is  found 
to  be  tractable. 


3.  Analysis  of  control  problem 

The  base  turbulent  flow  analyzed  in  this  problem  is  illustrated  in  figure  2.  This  flow 
has  been  carefully  studied  by  Kim,  Moin,  &  Moser  (1987);  the  present  simulations  yield 
essentially  the  same  second  order  statistics  as  the  results  presented  there  for  the  un¬ 
controlled  flow,  though  the  numerical  method  used  in  the  present  work  is  substantially 
different  in  order  to  better  facilitate  wall-normal  velocity  boundary  conditions. 

As  direct  numerical  simulations  of  turbulence  produce  a  tremendous  amount  of  data,  it 
is  important  to  analyze  relevant  statistics  of  these  flowfields  in  order  to  better  understand 
the  phenomena  taking  place  in  an  integrated  sense  and  how  these  integral  measures  of  the 
turbulence  are  modified  by  the  addition  of  control.  The  statistics  used  to  examine  the  tur¬ 
bulent  flowfields  in  the  present  work  are  the  mean  velocity  ui^  the  root-mean-square  ve¬ 
locity  fluctuations  Ui^rms^  fbe  Reynolds  stress  —uiU2^  the  total  stress  —uiU2  +  i'd^fdx2^ 
and  the  two-point  correlations  Rij{r)  ~  Ui{x)  Uj{x r)  and  their  Fourier  transform,  the 
cospectra  F?iy(k).  Note  that  the  overbar  (  )  implies  averages  in  the  homogeneous  direc¬ 
tions  xi  and  X2  and,  when  appropriate,  time.  The  statistics  are  functions  of  the  nonho- 
mogeneous  direction  X2 .  A  further  discussion  of  these  statistics  and  their  behaviour  in  a 
turbulent  channel  at  Rsr  =  180  may  be  found  in  Kim,  Moin,  &  Moser  (1987). 
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(a)  Three-dimensional  visualization.  For  clarity ,  discriminant  is  marked  only  in  the 
lower  half  of  the  domain.  Flow  is  from  left  to  right,  walls  are  shaded. 


rows )  are  marked  on  only  one  ninth  of  the  gridpoints  used  in  the  computation. 

Figure  2.  Turbulent  channel  flow  realization  at  Rcr  =  180,  no  control.  Regions  of  the  flow  with 
positive  discriminant  D  >  -Dthreshoid  are  shaded,  indicating  fluid  motion  which,  in  a  pointwise 
sense,  is  vortical  in  nature,  as  suggested  by  Blackburn,  Mansour,  &  Cantwell  (1996).  Small 
amounts  of  blowing  and  suction  will  be  applied  through  the  computational  equivalent  of  closely 
spaced  holes  drilled  in  the  walls  in  response  to  these  turbulent  motions  in  a  manner  which 
reduces  drag. 


The  nature  of  the  turbulent  motion  is  also  well  characterized  by  observing  the  fluid  at 
various  points  throughout  the  channel  in  a  reference  frame  which  moves  with  the  local 
velocity.  In  this  reference  frame,  the  point  under  consideration  is  a  critical  point,  as  the 
local  streamline  slope  is  indeterminate.  Thus,  a  critical  point  analysis  of  the  type  dis¬ 
cussed  by  Perry  &  Chong  (1987)  is  appropriate.  Chong,  Perry,  &  Cantwell  (1990)  and 
Blackburn,  Mansour,  &  Cantwell  (1996)  have  demonstrated  that  a  single  scalar  quantity 
D,  the  discriminant  of  the  velocity  gradient  tensor,  provides  a  useful  identification  of  re¬ 
gions  in  the  flow  which,  in  this  context,  are  “focus”  in  nature.  Such  focus  regions  roughly 
correspond  to  “vortex- type”  regions  in  a  turbulent  fiowfield,  though  this  description  is 
only  pointwise  in  nature. f 

The  velocity  gradient  tensor  discussed  in  this  work  is  defined  in  wall  units  Aij  = 
duf  jdx^ .  The  second  and  third  invariants  of  A  are  Q  =  {[tr(A)]^  ~  tr(^^)}  /2  and  = 
—  det(A).  The  discriminant  of  the  velocity  gradient  tensor  is  given  hy  D  (27/4)iZ^+Q^. 
Regions  with  D  >  0  are  characterized  by  a  velocity  gradient  tensor  with  one  real  and  two 
complex  eigenvalues  (and  thus  a  swirling,  vortex-type  motion  in  a  Lagrangian  reference 
frame),  whereas  regions  with  D  ^  0  are  characterized  by  three  real  eigenvalues.  For 
clarity,  the  visualizations  of  the  discriminant  presented  in  this  work  identify  only  regions 

f  Note  that  this  description  of  a  “vortex”  is  by  no  means  unique.  Robinson  (1991)  and 
Bernard,  Thomas,  &  Handler  (1993)  discuss  other  vortex  identification  techniques. 
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of  positive  discriminant  greater  than  a  small  threshold  value  D  >  /^threshold)  where 

■^threshold  —  10 


4.  Summary  of  optimal  control  theory 

4.1.  Cost  functional 

The  first  step  in  solving  an  optimal  control  problem  is  to  represent  the  control  problem 
of  interest  as  a  cost  functional,  to  be  minimized.  In  the  present  problem,  control  is 
to  be  applied  to  minimize  the  drag  averaged  over  a  section  of  wall  with  area  A  and  over 
the  time  period  (OjT]  using  the  least  amount  of  control  effort  possible.  A  relevant  cost 
functional  for  the  present  problem  is  thus 


j(#)  = 


dui 


dt  dS. 


The  first  term  in  the  integrand  is  a  measure  of  the  magnitude  of  the  control.  The  second 
term  is  a  measure  of  exactly  that  quantity  we  would  like  to  regulate:  in  this  case,  the 
drag.  These  quantities  are  integrated  over  the  wall  section  under  consideration,  of  area  A, 
and  over  the  time  period  under  consideration,  of  duration  T.  Finally,  they  are  weighted 
together  with  a  factor  £^/2,  which  represents  the  price  of  the  control.  This  quantity  is 
small  if  the  control  is  “cheap”  (which  reduces  the  significance  of  the  first  term) ,  and  large 
if  applying  control  is  “expensive”. 


4.2.  Gradient  of  cost  functional 

As  suggested  by  Abergel  and  Temam  (1990),  a  rigorous  procedure  may  be  developed  to 
determine  the  sensitivity  of  a  cost  functional  J  to  small  modifications  of  the  control  ^ 
for  nonlinear  problems  of  this  sort.  To  do  this,  consider  the  perturbation  to  the  cost 
functional  resulting  from  a  small  perturbation  to  the  control  ^  in  the  direction  .  (Note 
that  this  control  perturbation  direction  is  arbitrary  and  scaled  to  have  unit  norm.) 
Define  as  the  Frechet  differential  (Vainberg  1964)  of  the  cost  functional  such  that 


+  J(^) 

e-fO  e 


9J{^) 


^'dtdS. 


The  quantity  ff'  is  the  cost  functional  perturbation  due  to  a  control  perturbation 
scaled  by  the  inverse  of  the  control  perturbation  magnitude  e  in  the  limit  that  e  — )>  0.  The 
above  relation,  considered  for  arbitrary  also  defines  the  gradient  of  the  cost  functional 
J  with  respect  to  the  control  which  is  written 

In  the  present  approach,  the  cost  functional  perturbation  ff'  defined  above  is  expressed 
as  a  simple  linear  function  of  the  direction  of  the  control  perturbation  through  the 
solution  of  an  adjoint  problem.  By  the  above  formula,  such  a  representation  reveals  the 
gradient  direction  directly.  With  this  gradient  information,  the  control  ^ 

is  updated  on  (0,T]  in  the  direction  that,  at  least  locally  {i.e,  for  infinitesimal  control 
updates),  most  effectively  reduces  the  cost  functional.  The  finite  distance  the  control 
is  updated  in  this  direction  is  then  found  by  a  line  search  routine,  which  makes  this 
iteration  procedure  very  robust  even  when  controlling  nonlinear  phenomena.  The  flow 
resulting  from  the  modified  control  is  then  computed  according  to  the  (nonlinear)  Navier- 
Stokes  equation  (1.1),  the  sensitivity  of  this  new  flow  to  further  control  modification  is 
computed,  and  the  process  repeated.  Upon  convergence  of  this  iteration,  the  flow  is 
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advanced  over  the  interval  (OjTi],  where  Ti  ^  T,  and  an  iteration  for  the  optimal  control 
over  a  new  time  interval  (Ti ,  Ti  +  T]  begins  anew.f 

The  cost  functional  perturbation  resulting  from  a  control  perturbation  in  the  di¬ 
rection  is  given  by 


(4.1) 


where  is  the  Frechet  differential  of  wi,  as  defined  in  the  following  subsection.  Adjoint 
calculus  is  used  simply  to  re-express  the  integral  involving  as  a  linear  function  of 
Once  this  is  accomplished,  is  factored  out  of  the  integrands  and,  as  the  equation  holds 
for  arbitrary  an  expression  for  the  gradient  is  extracted. 


4.2.1.  Perturbation  field 

Define  [/'  as  the  Frechet  differential  of  U  such  that 

u’ .  Ii„  + 

e-+0  6 

The  quantity  U'  is  the  flow  perturbation  due  to  a  control  perturbation  scaled  by  the 
inverse  of  the  control  perturbation  magnitude  e  in  the  limit  that  e  — >  0.  The  equations 
governing  the  dependence  of  the  flow  perturbation  U'  on  the  direction  of  the  control 
perturbation  may  be  found  by  taking  the  Frechet  differential  of  the  state  equation 
(1.1)  itself.  The  result  is 


A/*'  C/'  =  0,  in  O 


(4.2a) 


with  boundary  conditions 


u-  =  Hi  on 


(4.2b) 


and  initial  conditions 


u[  =  0  at  /:  =  0, 

where  the  differential  Navier-Stokes  operator  fif'  is  given  by 


M'U'  = 


dt  ^  ^dxj  ^  ^dxj 
_1^ 
P  dxj 


,  1  \ 
'"'d^'^pdx 


(4.2c) 


(4.3) 


Note  that  the  operation  fif'  U'  is  linear  in  the  perturbation  field  t/',  though  the  operator 
J\f'  itself  is  a  function  of  the  solution  U  of  the  Navier-Stokes  problem.  Equation  (4.2) 
reflects  the  linear  dependence  of  the  perturbation  field  in  the  interior  of  the  domain 


t  Note  that  the  flow  may  be  advanced  over  a  time  interval  which  is  shorter  than  the  complete 
interval  over  which  the  optimization  was  performed.  The  rationale  for  such  an  approach  is  that 
the  controls  computed  near  the  end  of  each  optimization  interval  are  computed  without  regard 
to  the  (inevitable)  further  development  of  the  flow  beyond  the  end  of  the  optimization  interval. 
Thus,  the  controls  near  the  end  of  the  optimization  interval  are  not  as  effective  in  the  long  run  as 
those  controls  near  the  beginning  of  the  interval,  which  are  optimized  with  greater  “foresight” 
about  the  flow  development.  Thus,  the  controls  optimized  near  the  end  of  one  interval  may 
be  improved  upon  (in  terms  of  the  long-term  performance)  by  recalculation  in  the  following 
interval. 
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on  the  direction  of  the  control  perturbation  at  the  boundary.  However,  the  implicit 
linear  relationship  U'  =  [/'($')  given  by  these  equations  is  not  tractable  for  expressing 
J’'  in  a  form  from  which  may  be  deduced.  For  the  purpose  of  determining 

a  more  useful  relationship  with  which  we  may  express  ff'  in  the  desired  form,  we  now 
appeal  to  an  adjoint  identity. 


4.2.2.  Derivation  of  adjoint  identity 

This  subsection  derives  the  adjoint  of  the  linear  partial  differential  operator  Af' . 
Define  an  inner  product  over  the  domain  in  space-time  under  consideration  such  that 

([/',[/*)  =  f  [  U'  -U^dtdV, 

Jq  Jo 

and  consider  the  identity 


{Af'U',  U*)  =  {U',U'*U*)+  b. 


(4.4) 


Integration  by  parts  may  be  used  to  move  all  differential  operations  from  U'  on  the  left 
hand  side  of  the  equation  to  U*  on  the  right  hand  side,  resulting  in  the  definition  of  the 
adjoint  operator 


dur  /dui  duU  idp*\ 

dt  ^  V  dxj  dxi )  dx^j  p  dxi 

Idu] 

\  p  dxj  ) 


(4.5) 


where  the  operator  A/**  is  a  function  of  JJ ^  and  an  expression  for  5,  which  contains  all  the 
boundary  terms: 


The  identity  (4.4)  is  very  powerful;  in  fact,  simplification  of  this  identity  by  interior 
equations,  boundary  conditions,  and  initial  conditions  on  and  C/*  provides  an 

expression  which  recasts  in  a  form  from  which  may  be  deduced,  as 

illustrated  in  the  following  two  subsections. 


4.2.3.  Definition  of  adjoint  field 

Consider  an  adjoint  state  defined  (as  yet,  arbitrarily)  by 


A/**  =0,  in  (4.6a) 

with  boundary  conditions 

u*  =  Sii  on  (4.6b) 

and  initial  conditions 

w*  =  0  at  ^  T,  (4.6c) 

where  the  adjoint  operator  Af*  is  given  in  (4.5).  Note  that  the  adjoint  problem  (1.2), 
though  linear,  has  complexity  similar  to  that  of  the  Navier-Stokes  problem  (1.1).  Note 
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also  that  the  “initial”  conditions  in  (4.6c)  are  defined  at  i  =  T.  With  this  definition, 
the  adjoint  field  must  be  marched  backward  in  time  over  the  interval — due  to  the  sign 
of  the  time  derivative  and  viscous  terms  in  the  adjoint  operator  A/**  in  (4.5),  this  is  the 
natural  direction  for  this  time  march.  However,  as  the  operator  is  a  function  of  ff, 
computation  of  the  adjoint  field  U*  requires  storage  of  the  flow  field  U  ont  E  [0,  T],  which 
itself  must  be  computed  with  a  forward  march.  This  storage  issue  presents  a  numerical 
complication  which  precludes  solution  for  large  optimization  invervals  T.  However,  the 
problem  is  not  insurmountable  for  moderate  values  of  T. 

Equation  (1.2)  is,  as  yet,  simply  a  definition  of  an  adjoint  field.  The  motivation  for 
considering  an  adjoint  field  so  defined  is  revealed  in  the  following  subsection. 


4.2.4.  Identification  of  gradient 

The  identity  (4.4)  is  now  simplified  using  the  equations  defining  the  state  field  (1.1), 
the  perturbation  field  (4.2),  and  the  adjoint  field  (1.2).  Due  to  the  judicious  choice  of 
RHS  forcing  terms  in  the  equations  (1.2a)-(4.6c)  defining  the  adjoint  field,  the  identity 
reduces  to 


fi  dt  dS  — 
an 


p*  dtdS. 


Using  this  equation,  the  cost  functional  perturbation  fl'  in  (4.1)  may  be  conveniently 
rewritten  as 


f 

JuJo  \ 


dtdS==0. 


As  is  arbitrary,  this  implies  that 

(4.7) 

Thus,  the  desired  gradient  is  found  to  be  a  function  of  the  solution  of  the  adjoint  problem 
(1.2)  discussed  in  §4.2.3. 


4.3.  Gradient  update  to  control 

4.3.1.  Simple  gradient 

A  control  strategy  using  a  simple  gradient  (also  known  as  steepest  descent)  algorithm 
may  now  be  proposed  such  that 


_  ^k-l 


—  a 


(4.8) 


over  the  entire  time  interval  t  G  (0,T],  where  k  indicates  the  iteration  number  and  a  is 
a  parameter  of  descent  which  governs  how  large  an  update  is  made,  which  is  adjusted 
at  each  iteration  step  to  be  that  value  which  minimizes  ff.  This  algorithm  updates  $ 
at  each  iteration  in  the  direction  of  maximum  decrease  of  ff.  As  k  —>  oo^  the  algorithm 
should  converge  to  some  local  minimum  of  ff  over  the  domain  of  the  control  $  on  the 
time  interval  t  G  (0,T].  Note  that  convergence  to  a  global  minimum  will  not  in  general 
be  attained  by  such  a  scheme,  and  that,  as  time  proceeds,  ff  will  not  necessarily  decrease. 


4.3.2.  Conjugate  gradient 

The  simple  gradient  approach  described  above  is  straightforward,  but,  as  illustrated 
in  figure  3,  not  always  efficient.  Even  in  linear  problems,  for  cases  in  which  the  cost 
functional  has  a  long,  narrow  “valley”,  the  lack  of  a  momentum  term  from  iteration  to 
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iteration  tends  to  cause  the  simple  gradient  algorithm  to  bounce  from  one  side  of  the 
valley  to  the  other  without  turning  to  proceed  along  the  valley  floor.  The  conjugate 
gradient  approach  tends  to  improve  this  behaviour,  and  will  be  examined  further  in 
future  work. 

5.  Numerical  approach 

5.1.  Numerical  algorithm,  for  solution  of  flow  and  adjoint  equations 

The  control  formulations  derived  above  were  tested  in  direct  numerical  simulations  of 
fully  developed  turbulent  channel  flows  at  Rsr  =  180.  Fourier  transforms  are  used  to 
compute  spatial  derivatives  in  the  homogeneous  directions  with  3/2  dealiasing  on  the 
nonlinear  terms,  and  a  conservative  second-order  finite  difference  scheme  is  used  to  com¬ 
pute  spatial  derivatives  in  the  wall-normal  direction.  For  the  present  simulations,  the 
number  of  Fourier  modes  used  is  170  X  129  x  170  in  the  xi,  X2^  and  directions  respec¬ 
tively  {i.e.  256  X  129  X  256  dealiased  collocation  points),  and  the  size  of  the  computational 
domain  in  wall  units  is  L'^  =  2260,  Lj  =  360,  Lj  =  1130.  The  resulting  effective  grid 
resolution  in  the  streamwise  and  spanwise  directions  (on  collocation  points  determined 
without  the  extra  3/2  padding)  is  A.7;^  =  13,  Ar/;J  =  7.  Hyperbolic  tangent  streching  of 
the  grid  is  used  in  the  wall-normal  direction,  resulting  in  a  grid  spacing  of  Ax^  =0.3 
adjacent  to  the  wall  and  AxJ  =  5  in  the  center  of  the  channel.  Fine  grid  resolution  is 
required  near  the  wall  to  resolve  the  shear  layer;  the  mesh  is  fine  in  this  direction  even  up 
to  the  center  of  the  domain  because  the  second  order  difference  scheme  used  to  compute 
the  derivatives  in  this  direction  is  numerically  dissipative.  The  computational  grid  is 
staggered  in  the  wall-normal  direction  to  prevent  decoupling  of  the  even  and  odd  modes 
of  the  pressure. 

The  flow  is  advanced  in  time  using  an  explicit  low-storage  third-order  Runge-Kutta 
method  for  terms  involving  xi  and  x^  derivatives  and  an  implicit  Crank- Nicholson 
method  at  each  Runge-Kutta  substep  for  X2  derivative  terms.  A  temporal  discretiza¬ 
tion  implicit  in  the  X2  derivatives  is  necessary  to  mitigate  the  CFL  time  step  restriction 
when  control  is  applied,  as  the  control  fluid  at  the  wall  is  directed  in  the  X2  direction, 
which  is  precisely  the  region  and  direction  in  which  the  mesh  must  be  refined  most  to 
resolve  the  shear  layer. 

The  adjoint  solver  is  coded  with  a  method  analogous  to  that  of  the  flow  solver.  The 
flow  field  is  stored  every  5  time  steps  on  the  forward  sweep,  with  linear  interpolation  of 
these  stored  fields  used  on  the  backward  sweep  to  determine  the  operator  J\f* .  In  the 
optimal  calculations  presented  here,  we  chose  £  =  10“^  (control  power  is  taken  to  be 
“cheap”).  The  Polak-Ribiere  variant  of  the  conjugate  gradient  algorithm  was  used  for 
the  control  update,  with  a  computed  at  each  iteration  by  Brent’s  method,  a  robust  line 
minimization  algorithm  taken  from  Press  et  al  (1986). 


6.  Performance  of  controlled  systems 

At  the  time  this  document  went  to  press,  high  Reynolds  number  computational  results 
were  not  yet  available.  They  will  appear  in  the  version  of  this  paper  that  gets  submitted 
to  J.  Fluid  Mech.  Instead,  shown  here  are  some  of  the  visualizations  and  statistics  from 
the  baseline  case  (no  control)  which  was  run  in  order  validate  the  code  and  illustrate  the 
post  processing  analysis  to  be  performed. 

Figure  6  demonstrates  that  effective  performance,  including  about  50%  drag  reduction 
and  nearly  an  order  of  magnitude  reduction  in  the  turbulent  kinetic  energy,  have  been 
obtained  in  Rcr  =  100  turbulent  flow. 


iinant  in  unmanipulated  turbulent  channel!  now.  Areas  oi  positive  discriminant  JD  >  .^threshold  are  sJ 
indicating  fluid  motion  which,  in  a  pointwise  sense,  is  vortical  in  nature. 
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Figure  4.  Selected  instantaneous  statistics  of  the  unmanipulated  turbulent  flow. 
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(a)  History  of  drag. 


(b)  History  of  turbulent  kinetic  energy. 

Figure  5.  Performance  of  optimal  control  algorithm  at  Rcr  =  100: - ,  optimal 

drag  formulation  (§4);  - ,  optimal  energy  formulation; - ,  intuition  based  scheme 

$  =  — ^2(2/“*“  =  10);  . ,  no  control  (both  turbulent  and  laminar  drag  curves  are  shown). 
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PART  C. 

Robust  control  of  turbulence 

The  framework  for  applying  optimal  and  robust  control  theories  to  linear  problems 
is  first  reviewed  in  a  notation  fairly  consistent  with  that  used  by  Doyle  et  al  (1989). 
The  discussion  is  made  strictly  in  the  time  domain,  not  the  frequency  domain  often 
used  to  discuss  these  approaches,  to  facilitate  extension  of  the  approaches  to  nonlinear 
problems  in  which  frequency  domain  techniques  are  of  limited  usefulness.  The  resulting 
development  is  fairly  straightforward  and  does  not  assume  the  reader  is  accustomed  to 
the  language  of  control  theory. 

Nonlinear  optimal  control  problems  are  then  reviewed  in  a  similar  notation— the  ap¬ 
proach  described  here  is  that  taken  by  Bewley,  Moin,  and  Temam  (1997)  for  the  optimal 
control  of  wall-bounded  turbulent  flows.  Finally,  the  concepts  of  robust  control  are  ex¬ 
tended  to  nonlinear  problems,  such  as  the  control  of  turbulence,  in  a  consistent  manner. 


1.  Outline  of  linear  regulation  problems 

1.1.  Optimal  regulation  of  linear  problems 

1.1.1.  State  equation 

Consider  a  state  vector  u  which  is  a  function  of  some  feedback  control  vector  ^  such 
that  it  obeys  the  linear  evolution  equation 


u=  Au+  82^ 


(l.la) 


with  given  initial  conditions 


u  =  u(0)  at  ^  =  0.  (1-lb) 

The  matrices  A  and  B2  may  be  functions  of  time  but  do  not  themselves  depend  on  the 
state  u  or  the  control 


1.1.2.  Cost  function 

The  object  of  applying  control  in  the  present  problem  is  to  regulate  some  measure  of 
the  state  to  zero  quickly  without  applying  excessive  amounts  of  control.  Mathematically, 
this  objective  is  expressed  as  the  minimization  of  a  cost  functional  which  balances  a 
measure  of  the  state  u  with  a  measure  of  the  control  $  applied.  We  will  use  the  norm 
symbol  to  denote  these  measures,  which  may  be  defined  appropriately  for  particular 
problems  of  interest: 

Note  that  the  two  terms  are  weighted  together  with  a  factor  which  accounts  for  the 
“price”  of  the  control;  this  factor  is  large  if  applying  the  control  is  “expensive”,  which 
emphasizes  the  importance  of  the  latter  term  in  this  equation  and  generally  results  in 
a  modest  control  effort,  and  small  if  applying  the  control  is  “cheap”,  which  results  in  a 
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larger  control  effort.  These  terms  are  then  averaged  over  some  optimization  interval  of 
consideration  [0,T].  In  matrix  form,  J2  is  expressed  as 

=  +  dt 

with  C\  defined  appropriately  based  on  the  definition  of  the  norm  of  the  state  ||u||  and 
the  star  (*)  denoting  the  conjugate  transpose.  By  appropriate  scaling  of  the  vector  ^ 
and  the  matrix  B25  the  norm  of  the  control  ||#||  is  taken  simply  as  the  Euclidian  norm 
without  loss  of  generality. 

A  technique  to  design  a  feedback  control  relationship  of  the  form  ^  =  K2M  which 
minimizes  the  cost  function  J2  is  now  briefly  outlined. 

1.1.3.  Adjoint  equation 

Define  an  adjoint  state  (as  yet,  arbitrarily)  by  the  relation 

~A  =  A*A  +  Ci*C'iU  (1.2a) 


with  initial  conditions 


A  =  0  att  =  T.  (1.2b) 

Note  that  the  “initial”  conditions  (1.2b)  are  defined  at  ^  =  T,  so  to  determine  the 
adjoint  on  the  interval  [0,T),  the  evolution  equation  (1.2a)  must  be  marched  backv)ards 
from  T  — 0. 


1.1.4.  Gradient  of  cost  function  J2  with  respect  to  control  ^ 

It  may  be  shown  that  the  gradient  of  the  cost  with  respect  to  the  control  is  a  simple 
function  of  the  adjoint  state  defined  by  (1.2): 

^  =  +  (1.3) 


1.1.5.  Solution  of  control  problem 
By  (1.3),  the  most  suitable  control  which  results  in 

=  0  (minimum)  (1-4) 

as  a  function  of  the  adjoint  state  is  given  simply  by 


(1.5) 


Combining  the  state  equation  (1.1a),  the  adjoint  equation  (1.2a),  and  the  control  given 
by  (1.5)  into  a  combined  matrix  form  gives 


(1.6) 


Now  prescribe  a  relationship  between  any  state  vector  u  =  u(^)  and  the  corresponding 
adjoint  X  =  X{t)  such  that 


A  =  X2U, 


(1.7) 
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where  X2  =  X2{t).  Inserting  this  expression  into  (1.6)  to  eliminate  A  and  combining  the 
top  and  bottom  rows  to  eliminate  li  leads  to  the  expression 

(^-X2  =  A*  X2  +  X2  A- X2^B2B;X2  +  C;Ci)u 

As  this  expression  is  valid  for  any  state  vector  u,  we  arrive  at  a  Riccati  equation  for  the 
matrix  X2{t): 


-X2  =  A*X2+X2A-X2^B2  B;  X2  +  Cl  Cl  (1.8a) 

with  initial  conditions,  due  to  (1.2b)  and  (1.7),  given  by 

X2  =  0  at  f  =  T.  (1.8b) 

Combining  (1.5)  and  (1.7),  the  optimal  control  ^  as  a  function  of  the  state  u  is  given  by 
the  state  feedback  relationship 


1  , 

^  K2U 

where 

K2  =  -pBlX2 

(1.9) 


where  X2(t)  is  the  solution  of  (1.8)  and  thus  K2  =  X2(t). 


1.1.6.  Infinite  time  horizon  for  tim>e  invariant  problem.s 

If  the  matrices  A,  B2,  and  Ci  are  time  invariant,  then  in  the  limit  of  large  optimization 
intervals  T  — )•  00  the  matrix  X2{f)  defined  by  (1.8)  approaches  a  steady  state  value  in 
the  march  from  the  initial  conditions  defined  t  =  T  back  towards  t  =  0.  This  steady 
state  value  may  be  found  by  setting  X2  =  0  in  (1.8a),  which  leads  to 


o  =  a*X2-¥X2A-X2^B2b;x2  +  ci  Ci 


(1.10) 


The  optimum  feedback  relationship  given  by  (1.9)  in  this  limit  is  thus  time  invariant  and 
a  function  of  the  solution  to  (1.10),  referred  to  as  an  algebraic  Riccati  equation  (ARE). 
Solution  methods  for  equations  of  this  type  are  well  developed  (Laub  1991). 


1.2.  Simple  interpretation  of  the  adjoint  fi.eld 

In  the  preceding  discussion,  the  determination  of  optimal  feedback  control  relationship 
$  =  K2  u  in  (1.9)  was  closely  linked  to  the  definition  of  an  adjoint  A  in  (1.2).  However, 
the  definition  of  A  was  made  arbitrarily  in  (1.2),  and  subsequently  justified  only  mathe¬ 
matically  in  (1.3)  as  being  that  field  which  is  required  to  express  the  gradient  of  the  cost 
function  with  respect  to  the  control  in  a  simple  manner. 
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In  the  case  that  the  control  $  enters  the  state  equation  (1.1)  through  the  identity 
matrix,  B2  ^  ^  simple  interpretation  of  the  adjoint  is  now  clear.  In  this  case,  the 

expression  for  the  gradient  (1.3)  reduces  to 

^  =  £2$  + A.  (1.11) 

Thus,  the  gradient  consists  of  two  terms.  The  first  term  simply  accounts  for  the  term 
in  the  cost  function  J2  which  measures  the  magnitude  of  the  control;  in  the  absense  of 
other  term  in  the  cost  function,  this  term  would  drive  the  control  to  zero  when  J2  is 
minimized. 

The  second  term  in  (1.11)  accounts  for  the  term  in  the  cost  function  J2  which  measures 
the  state  u  itself.  Thus,  one  interpretation  of  the  adjoint  A  is  simply  that:  The  adjoint, 
'when  properly  defined,  is  a  measure  of  the  sensitivity  of  the  term  of  the  cost  function 
which  measures  the  state  u  to  additional  RHS  forcing  of  the  state  equation.  Note  that 
there  are  exactly  as  many  components  of  the  adjoint  A  as  there  are  components  of  the 
state  equation  (1.1a). 


1.3.  Robust  regulation  of  linear  problems 

1.3.1.  State  equation 

Consider  the  linear  state  equation  of  (1.1)  with  additional  forcing  due  to  an  external 
disturbance  x 


u  =  Au  -\-  Bi  X  B2  ^ 


(1.12a) 


with  given  initial  conditions 

u  =  u(0)  at  i  =  0.  (1.12b) 

The  matrix  Bi  may  be  a  function  of  time  but  does  not  itself  depend  on  the  state  u  or 
the  control 


1.3.2.  Cost  function 

The  object  of  applying  control  in  the  robust  problem  is  identical  to  the  optimal  prob¬ 
lem,  except  we  now  play  the  “devil’s  advocate”  and  seek  to  find  the  best  control  in  the 
presence  of  a  “small”  component  of  exactly  that  disturbance  %  which  is  maximally  ag- 
grevating  to  the  control  objective.  To  represent  this  concept  mathematically,  we  append 
to  the  cost  function  discussed  in  the  previous  section  a  term  which  accounts  for  the 
magnitude  of  the  disturbance  used  to  aggrevate  the  system 

Note  that  the  sign  of  the  term  which  is  used  to  account  for  the  disturbance  is  opposite 
to  the  sign  used  to  account  for  the  control;  this  is  because  we  minimize  with  respect 
the  control  $  while  simultaneously  we  m>aximize  with  respect  to  the  disturbance  %.  The 
term  involving  —7^  ||x|P  limits  the  magnitude  of  the  disturbance  in  the  maximization 
with  respect  to  x  the  term  involving  ||^||^  limits  the  magnitude  of  the  control  in 
the  minimization  with  respect  to  In  matrix  form,  Joo  is  expressed  as 

Joo  =  C{CiU  +  e  ^  x)  dt 
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By  appropriate  scaling  of  the  vector  x  fhe  matrix  Bi ,  the  norm  of  the  disturbance 
llxll  is  taken  simply  as  the  Euclidian  norm  without  loss  of  generality. 

A  technique  to  design  a  feedback  control  relationship  of  the  form  $  =  ATcjo  u  which 
minimizes  the  cost  function  j7oo  hi  the  presence  of  a  small  component  of  the  worst  external 
disturbance  x  forcing  the  state  equation  (1.12)  is  now  briefly  outlined.  By  designing  a 
feedback  control  rule  effective  for  a  state  disturbed  in  this  manner,  the  control  rule  which 
is  found  is  effective  in  the  presence  of  small  disturbances  of  any  type,  and  has  nearly  the 
same  nominal  performance  {i.e.  performance  on  the  undisturbed  system)  as  the  optimal 
controller  determined  in  the  previous  section. 

1.3.3.  Adjoint  equation 

Define  an  adjoint  state  as  for  the  optimal  control  case  by  the  relation 

-A  =  A*A  +  C*CiU  (1.13a) 

with  initial  conditions 

A  =  0  att  =  T.  (1.13b) 

1.3.4.  Gradients  of  cost  function  Joo  '^dth  respect  to  control  ^  and  distm'bance  x 

In  a  manner  identical  to  the  derivation  leading  to  (1.3),  the  gradient  of  the  cost  with 
respect  to  the  control  ^  and  the  disturbance  x.  fhis  problem  are  simple  functions  of 
the  adjoint  state  defined  by  (1.13): 

+  and  ^  =  +  (1.14) 


1.3.5.  Solution  of  control  problem 

By  (1.14),  the  most  suitable  control  and  disturbance  which  result  in 


c\  f  '  '  \  A  ^ 

— - —  =  0  (minimum)  and  — ^ —  —  0  (maximum) 

9^  ^  ^  !]^X 


(1.15) 


are  given  simply  by 

^  =  and  (1.16) 

Combining  the  state  equation  (1.12a)  and  the  adjoint  equation  (1.13a)  with  the  control 
and  disturbance  given  by  (1.16)  into  a  combined  matrix  form  gives 


~ 

1  .1  J 

u 

\ 

— 

A 

^ 

/\ 

_-ci  cr 

-A* 

u 

A 


(1.17) 


Now  prescribe  a  relationship  between  any  state  vector  u  and  the  corresponding  adjoint 
A  such  that 


A  =  XooU 


(1.18) 


Inserting  this  expression  into  (1.17)  to  eliminate  A  and  combining  the  top  and  bottom 
rows  to  eliminate  li  leads  to  the  expression 


-Xoc  =  A*  Xoc  +  Xoc  A  +  Xo 


Bi  B 


Xoo  +  Cl  Cl 


u 
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As  this  expression  is  valid  for  any  state  vector  u,  we  arrive  at  a  Riccati  equation  for  the 
matrix  Xoo{^)- 

-x^  =  A*x^  +  x^A  +  Xo^  (^Bib;-^B2b;]  x^  +  c;ci 

W  P  )  (1.19a) 

with  initial  conditions,  due  to  (1.13b)  and  (1.18),  given  by 

Xoc  =  0  ati  =  T.  (1.19b) 

Combining  (1.16)  and  (1.18),  a  robust  control  ^  which  is  effective  even  in  the  presence 
of  a  small  component  of  the  worst  case  disturbance  x  is  given  by  the  state  feedback 
relationship 


1  , 

where 

Koc  =  -^b;Xoo 

where  Xoo  is  the  solution  of  (1.19)  and  thus  K^o  =  R^ooi^)- 


1.3.6.  Infinite  time  horizon  for  time  invariant  problems 

If  the  matrices  A,  5i,  jB2,  and  Ci  are  time  invariant,  then  in  the  limit  that  the 
optimization  interval  T  — >  oo  the  matrix  Xoo(f)  defined  by  (1.19)  approaches  a  steady 
state  value  in  the  march  from  the  initial  conditions  defined  at  i  =  T  back  towards  t  =  0, 
and  is  given  by  the  solution  to 


0  =  A*  Xoc  +  XooA  +  Xoo 


Xoo  +  Cl*  Cl 


(1.21) 


The  robust  feedback  relationship  given  by  (1.20)  in  this  limit  is  thus  time  invariant  and 
a  function  of  the  solution  to  (1.21). 


2.  Outline  of  nonlinear  regulation  problems 

2.1.  Optimal  regulation  of  nonlinear  problems 

2.1.1.  Stale  equation 

Consider  a  state  vector  u  which  is  a  function  of  some  feedback  control  vector  ^  such 
that  it  obeys  the  nonlinear  evolution  equation 

(2.1a) 


with  given  initial  conditions 

u=u(0)  at  t  =  0.  (2.1b) 

The  nonlinear  functions  A(u)  and  B2{^)  may  themselves  be  functions  of  time. 

2.1.2.  Cost  function 

The  object  of  applying  control  in  the  present  case  is  identical  to  the  optimal  linear  reg¬ 
ulation  problem  described  in  §1.1,2.  Mathematically,  it  is  expressed  as  the  minimization 
of 

T 


U  =  4(u)  +  B2($) 


(2.2) 
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A  technique  to  determine  the  control  $  on  the  interval  (0,  T]  which  (locally)  minimizes 
the  cost  function  J2  for  the  nonlinear  state  equation  (2.1)  is  now  briefly  outlined. 


2.1.3.  Perturbation  equation 

Consider  the  linear  problem  of  a  small  perturbation  (^',u')  to  some  reference  solution 
($,u)  of  the  system  given  by  (2.1).  It  is  easily  shown  that  such  a  perturbation  must 
obey  a  linear  evolution  equation  of  the  form 


u'  =  A  u'  +  J32 


(2.3a) 


with  initial  conditions 


u'  =  0  at  i  =  0.  (2.3b) 

The  matrices  A  and  B2  are  functions  of  time  and  depend  explicitly  on  the  reference 
condition  (^jU). 

2.1.4.  Adjoint  equation 

An  adjoint  system  is  defined  based  on  the  A  matrix  of  the  perturbation  problem  (2.3) 
such  that 

=  A*A  +  Ci*Ciu  (2.4a) 

with  initial  conditions 

A  0  at  ^  T.  (2.4b) 

2.1.5.  Gradient  of  cost  function  J2  with  respect  to  control  $ 

As  in  the  linear  case,  the  gradient  of  the  cost  with  respect  to  the  control  is  a  simple 
function  of  the  adjoint  state  defined  by  (2.4): 

^  =  £2  $  +  b;  a.  (2.5) 

^$2  V  / 

2.1.6.  Solution  of  control  problem 

The  most  suitable  control  on  (0,  T]  which  results  in 

^^  =  0  (minimum)  (2.6) 

may  not  be  found  simply  by  setting  the  gradient  in  (2.5)  equal  to  zero,  as  this 

gradient  information  is  accurate  only  in  a  small  neighborhood  of  the  reference  solution 
upon  which  the  matrices  A  and  B2  were  based. f  Instead,  a  more  stable  iterative  approach 
is  used  based  on  the  gradient  vector: 

(2-^) 

where  k  indicates  the  iteration  index.  Thus,  the  condition  (2.6)  is  approached  iteratively 
according  to  the  following  procedure: 

t  One  may  propose  a  Newton-Raphson  technique  to  determine  the  control,  setting  the  local 
expression  for  in  (2.5)  equal  to  zero  to  determine  a  new  control,  determining  a  new 

reference  condition  from  (2.1),  examining  the  new  perturbation  problem  to  determine  a  new 
expression  for  and  iterating  until  convergence.  However,  such  a  technique  is  not 

recommended,  as  there  is  no  way  to  insure  that  the  initial  reference  condition  is  sufficiently 
close  to  a  minimum  to  guarantee  convergence  of  this  approach. 
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(а)  Initialize  control  ^  on  (0,T]  to  ^  =  0. 

(б)  Determine  state  u  on  (0,T]  from  state  equation  (2.1). 

(c)  Determine  adjoint  A  on  [0,T)  from  adjoint  equation  (2.4). 

(d)  Determine  local  expression  for  gradient  from  (2.5). 

(e)  Test  various  different  values  for  a  in  (2.7),  computing  the  resulting  state  u  from 
(2.1)  and  the  resulting  cost  J2  from  (2.2),  and  determine  via  a  line  minimization  algo¬ 
rithm  that  value  of  a  which  results  in  the  smallest  J2- 

(/)  Update  control  ^  on  (0,T]  via  (2.7)  using  best  value  of  a  determined  in  step  e. 

{g)  Repeat  from  step  b  until  convergence. 


2.2.  Robust  regulation  of  nonlinear  problems 

2.2.1.  State  equation 

Consider  the  nonlinear  state  equation  of  (2.1)  with  additional  forcing  due  to  an  external 
disturbance  x 


u=A{u)  +  Bi{x)  +  B2{^) 


(2.8a) 


with  given  initial  conditions 

u  =  u(0)  at  ^  0.  (2.8b) 

The  nonlinear  function  Bi  (x)  itself  be  a  function  of  time. 

2.2.2.  Cost  function 

The  object  of  applying  control  in  the  present  case  is  identical  to  the  robust  linear  reg¬ 
ulation  problem  described  in  §1.2.2.  Mathematically,  it  is  expressed  as  the  minimization 
of  a  cost  function  J^oo  with  respect  to  the  control  ^  while  simultaneously  maximizing  ffoo 
with  respect  to  the  disturbance  %,  where 


A  technique  to  determine  the  control  ^  on  the  interval  (0,T]  which  (locally)  mini¬ 
mizes  the  cost  function  ffoo  in  the  presence  of  a  small  component  of  the  worst  external 
disturbance  x  forcing  the  state  equation  (2.8)  is  now  briefly  outlined. 


2.2.3.  Perturbation  equation 

Consider  the  linear  problem  of  a  small  perturbation  (^',x^u')  to  some  reference  so¬ 
lution  (•i^jXjU)  of  the  system  given  by  (2.8).  It  is  easily  shown  that  such  a  perturbation 
must  obey  a  linear  evolution  equation  of  the  form 


u'  =  Au  Bi  x'  +  B2 


(2.9a) 


with  initial  conditions 


u'  =  0  at  f  =  0.  (2.9b) 

The  matrices  A,  ,  and  B2  are  functions  of  time  and  depend  explicitly  on  the  reference 
condition  ($,X)U). 
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2.2.4.  Adjoint  equation 

An  adjoint  system  is  defined  based  on  the  A  matrix  of  the  perturbation  problem  (2.9) 
such  that 

~X  =  A^  X  +  C^Ciu  (2.10a) 

with  initial  conditions 

A  =  0  at  t  =  T.  (2.10b) 

2.2.5.  Gradients  of  cost  function  ffoo  with  respect  to  control  $  and.  disturbance  x 

As  in  the  linear  case,  the  gradients  of  the  cost  with  respect  to  the  control  ^  and  the 
disturbance  x.  ^.re  simple  functions  of  the  adjoint  state  defined  by  (2.10): 

+  and  ^  =  +  (2.11) 


2.2.6.  Solution  of  control  problem. 

The  most  suitable  control  and  disturbance  which  result  in 

r\  (  '  '  \  J  ^ffoo  ^  /'o  1  o\ 

— - —  =  0  (minimum)  and  — - —  =  0  (maximum)  (2.12) 

may  not,  strictly  speaking,  be  found  simply  by  setting  the  gradients  and 

^Jool^X  in  (2.11)  equal  to  zero,  as  this  gradient  information  is  accurate  only  in  a  small 
neighborhood  of  the  reference  solution  upon  which  the  matrices  A,  .Bi,  and  B2  were 
based.  Instead,  an  iterative  approach  is  used  based  on  the  gradient  vectors: 

and  =  x"  + (2.13) 

^x 

where  k  indicates  the  iteration  index.  The  iteration  procedure  followed  is  analagous  to 
that  described  in  §2.1.6;  in  the  present  case,  a  value  of  a  is  chosen  to  reduce  Jc)o  while 
simulataneously  a  value  of  /5  is  chosen  to  increase  Joe-  The  min/max  problem  is  solved 
when  the  conditions  given  in  (2.12)  are  approached. 

2.2.7.  Approxim.ate  solution  for  systems  of  very  large  d.im.ension 

The  min/max  problem  described  by  (2.13)  is  unfeasible  when  the  state  equation  (2.8) 
is  a  model  of  turbulent  channel  flow,  as  the  state  u  upon  which  the  disturbance  acts  in 
this  case,  and  therefore  any  general  representation  of  the  disturbance  x.  itself,  has  a  very 
large  dimension  (C?(10^)  at  Rcr  —  180).  Thus,  instead  of  forcing  the  state  equation  with 
a  disturbance  x  determined  by  the  iterative  approach  given  in  (2.13),  which  is  guaranteed 
to  be  stable  but  would  present  excessive  computational  storage  requirements,  we  settle 
on  a  simpler,  though  possibly  unstable,  approach  for  the  determination  of  x.f  To  do 
this,  we  choose  the  noise  by  setting  ^Jooj^X.  in  (2.11)  equal  to  zero  to  determine  the 
the  disturbance  %.  Taking  the  matrix  Bi  as  simply  the  identity  matrix,  the  disturbance 
determined  in  this  fashion  is  proportional  to  the  adjoint  field  itself 

(2.14) 

f  Note  that  we  still  determine  the  control  $  via  the  stable  iterative  approach  given  in  (2.13). 
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For  sufficiently  large  7  sufficiently  small  level  of  noise),  this  is  an  accurate  ap¬ 

proximation  of  the  global  maximum  S^Joo/^X  =  0,  and  thus  results  in  an  accurate 
approximation  of  the  “worst  case”  noise.  For  smaller  values  of  7  larger  noise  lev¬ 
els),  this  approach  can  not  even  be  guaranteed  to  be  stable.  Trial  and  error  will  indicate 
for  what  values  of  7  this  approach  converges. 
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PART  D. 

Optimization  of  practical  feedback  rules 
for  turbulence  control 

A  new  method  based  on  control  theory  for  optimizing  feedback  control  rules  with  the 
objective  of  reducing  drag  in  wall-bounded  turbulent  flows  is  presented.  Both  linear  and 
nonlinear  control  rules  (of  the  type  commonly  used  in  neural  networks)  are  considered. 
These  control  rules  relate  wall  measurements  of  skin  friction  and  pressure  to  the  control, 
which  is  applied  as  a  continuous  distribution  of  wall-normal  boundary  velocity  with  zero 
net  transpiration.  Though  the  optimization  technique  itself  requires  complete  informa¬ 
tion  about  the  flow,  and  thus  can  only  be  performed  computationally,  it  is  intended  that 
the  resulting  optimized  rules  be  scaled  appropriately  and  used  in  physical  boundary  layer 
control  implementations. 

Using  optimal  control  theory,  the  sensitivity  of  some  representative  cost  functional 
to  small  modifications  in  the  coefficients  of  a  feedback  control  rule  are  found  via  the 
solution  of  an  adjoint  problem.  With  this  sensitivity  field,  the  coefficients  are  iteratively 
updated  with  a  gradient  algorithm  until  the  cost  functional  is  minimized.  Given  that  this 
optimization  is  performed  in  a  representative  situation,  the  coefficients  then  may  be  fixed 
and  the  control  rule  effectively  used  in  other  flows  with  similar  configurations,  requiring 
only  information  about  the  flow  which  can  be  obtained  with  flush-mounted  sensors  on 
the  wall. 


1.  Introduction 

Optimal  control  theory  applied  to  turbulence  provides  a  rigorous  framework  to  deter¬ 
mine  the  gradient  of  a  cost  functional  (which  represents  a  physical  problem  of  interest) 
with  respect  to  small  modifications  of  the  control  forcing  (Bewley,  Temam,  and  Moin, 
1997).  With  such  information,  combined  with  a  gradient  algorithm  to  update  the  con¬ 
trol,  very  effective  control  distributions  may  be  determined.  Numerical  simulations  of 
this  approach  in  a  low  Reynolds  number  turbulent  channel  flow  obtained  a  50%  drag 
reduction  and  an  order  of  magnitude  turbulent  kinetic  energy  reduction  with  small  levels 
of  boundary  velocity  control.  Important  drawbacks  of  this  approach,  however,  are  1)  it 
requires  complete  information  about  the  turbulent  fluctuations  in  the  near-wall  region, 
and  2)  it  is  extremely  computationally  expensive.  Thus,  it  is  impossible  to  apply  the 
optimal  control  approach  directly  in  an  experimental  setting. 

In  order  to  arrive  at  a  practical  scheme,  a  method  was  sought  to  optimize  control 
rules  which  1)  require  only  flow  information  obtainable  with  wall-mounted  sensors,  and 
2)  are  computationally  inexpensive  enough  to  apply  in  real  time.  Possible  approaches 
for  this  purpose  can  be  divided  into  two  broad  categories:  state  trajectory  approaches^ 
which  attempt  to  drive  some  description  of  the  turbulent  state  (or  a  portion  thereof)  in 
a  desired  manner,  and  direct  approaches^  which  bypass  any  description  of  the  turbulent 
state  per  se,  but  simply  seek  a  control  rule  which  achieves  a  desired  effect,  such  as  the 
reduction  of  drag. 

As  an  example  of  one  state  trajectory  approach,  an  adaptive  inverse  technique  has 
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been  applied  to  a  low  Reynolds  number  turbulent  channel  flow,  providing  approximately 
18%  drag  reduction  (Kim,  1996).  This  approach  first  develops  an  approximate  “inverse” 
model  between  measurable  flow  quantities  (as  input)  and  the  control  forcing  (as  output) 
with  an  adaptive  technique.  Each  iteration  of  the  adaptation  consists  of  three  steps:  1) 
computing  the  error  of  the  model  output  with  respect  to  the  desired  model  output  (the 
actual  control  forcing  used),  2)  determining  the  influence  of  the  weights  in  the  model  on 
this  error,  then  3)  updating  all  the  weights  in  the  model  a  small  amount  in  a  manner  that 
reduces  the  error.  In  neural  networks,  this  is  commonly  referred  to  as  “back- propagation” 
of  the  error.  Once  this  approximate  inverse  model  between  the  flow  measurements  and 
the  control  converges,  the  inverse  model  is  used  to  compute  a  control  which  will  drive  the 
flow  measurements  to  some  desired  state.  In  the  case  of  Kim  (1996),  the  desired  state  is 
chosen  to  be  a  state  with  reduced  spanwise  drag  fluctuations. 

Drawbacks  of  the  adaptive  inverse  approach  are  1)  an  ad  hoc  desired  state  must  be 
chosen,  2)  a  random  “dither”  signal  needs  to  be  applied  to  the  control  in  order  for  the 
inverse  model  to  have  “persistently  exciting”  data  from  which  to  learn,  which  reduces 
the  performance  of  the  controller,  and  3)  it  is  possible  that  even  at  statistical  steady 
state,  due  to  the  nonlinear  nature  of  the  Navier-Stokes  equation,  the  weights  in  the 
inverse  model  may  need  to  continually  adapt  in  order  to  represent  a  temporally  evolving 
relationship  between  the  flow  measurements  and  the  control.  Thus,  if  the  training  of  the 
inverse  model  does  not  converge  fast  enough,  it  will  not  have  time  to  keep  up  with  the 
temporal  evolution  of  the  flow  (for  instance,  the  movement  of  the  near- wall  turbulent 
coherent  structures),  and  may  not  develop  an  accurate  model  between  flow  measurements 
and  the  control  which  produces  them. 

Other  state  trajectory  approaches  attempt  to  control  a  more  complete  description  of 
the  turbulence  using  a  low-dimensional  (10-20  mode)  representation  of  the  near-wall 
coherent  structures  (Coller  et  al.^  1994).  In  this  approach,  the  orbit  of  the  near  wall 
structures  in  this  representation  is  partially  stabilized,  resulting  in  a  reduced  “bursting” 
frequency  and,  presumably,  reduced  drag,  Coller  et  aL  (1994)  showed  that  the  frequency 
of  bursting  events  could  be  reduced  in  their  model  equations,  but  did  not  demonstrate 
how  effective  such  an  approach  would  be  at  reducing  drag  when  applied  to  a  fully  tur¬ 
bulent  flow. 

Drawbacks  of  this  low- dimensional  representation  approach  include  1)  an  accurate 
estimation  of  this  low- dimensional  state  needs  to  be  made  from  the  measurements  at 
the  wall,  and  2)  a  desired  ad  hoc  state  trajectory  must  be  chosen,  which  can  only  be 
selected  well  if  one  has  a  detailed  understanding  of  the  cause/effect  relationship  of  the 
drag-producing  phenomena  in  the  near-wall  region,  which  is  still  under  debate. 

Direct  approaches  may  be  proposed  which  bypass  estimation  and  control  of  the  state 
trajectory  altogether.  In  such  approaches,  one  simply  represents  the  control  objective 
mathematically  as  a  cost  functional,  then  attempts  to  find  a  control  rule  which  minimizes 
this  functional. 

The  simplest  direct  approach  is  an  adaptive  “reinforcement  learning”  approach.  In 
such  an  approach,  the  weights  of  a  control  rule  are  initialized  randomly  and  the  control 
rule  applied  to  the  flow.  Every  time  a  “good”  result  is  seen  {e.g.,  the  drag  is  reduced), 
the  weights  contributing  most  to  the  control  at  that  instant  are  increased,  and  every  time 
a  “bad”  result  is  seen,  the  corresponding  responsible  weights  are  decreased. 

The  main  drawback  of  this  approach,  however,  is  that  this  reward/punish  training 
algorithm  is  not  very  reliable,  especially  for  complicated  nonlinear  systems,  and  thus  the 
scheme  may  not  converge  at  all. 

Thus,  we  arrive  at  the  motivation  for  the  current  work,  in  which  we  derive  a  rigorous 
algorithm  to  efficiently  optimize  a  direct  control  scheme,  with  the  goal  in  mind  simply 
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of  reducing  some  integral  measure  of  the  control  objective  without  the  prescription  of 
a  desired  state  trajectory.  This  approach,  based  on  computation  of  the  gradient  of  a 
cost  functional  with  respect  to  modification  of  the  weights  in  the  control  rule,  will  be 
outlined  in  the  following  sections.  Numerical  simulations  that  implement  this  technique 
are  underway. 


2.  Problem  Statement 


Our  goal  is  to  determine  a  relationship  which  takes  as  input  the  measurable  flow 
quantities  and  produces  as  output  a  control  $  (the  normal  component  of  velocity  at  the 
wall)  which  effectively  controls  the  flow  system.  The  measurable  flow  quantities  are  taken 
to  be  the  wall  values  of  the  streamwise  drag  fjiduijdn^  the  pressure  p,  and  the  spanwise 
drag  iJidu^ldrij  where  Ui  are  the  components  of  the  velocity,  /i  is  the  viscosity,  and  rii 
is  a  unit  normal  on  each  wall  facing  into  the  flow.  Hence,  one  of  the  simplest  (linear) 
control  rules  which  may  be  proposed  for  this  purpose  is 


dui  du3 

—  wi*  fj,— - h  W2 

on  on 


(2.1a) 


The  task  at  hand  is  simply  to  optimize  the  weighting  functions  w^,  such  that  the  control 
$  determined  by  linear  relationship  (2.1a)  effectively  controls  the  flow  system.  This 
configuration  is  illustrated  graphically  in  Figure  1.  Note  that  the  w^.  are  convolution 
functions,  where  the  convolution  is  defined  such  that,  for  example, 

j  Wi{x)  -  x)  dr,  (2.1b) 


where  T  is  the  portion  of  the  (2D)  boundary  of  the  (3D)  flow  domain  ^  over  which 
measurements  are  made  (and  also,  we  assume,  the  portion  of  the  boundary  over  which 
control  is  applied),  and  x  G  f  is  the  variable  of  integration.  By  optimizing  the  convolution 
functions  ty«;(x),  we  take  into  account  “nearby”  flow  measurements  (in  the  direction  x) 
from  a  specific  actuator  location  {x).  In  fact,  the  extent  to  which  these  convolution 
functions  are  nonzero  when  converged  will  indicate  how  far  in  each  direction  from  a 
specific  actuator  flow  measurements  are  relevant  when  computing  an  effective  control. 
Additionally,  we  will  constrain  the  spatial  average  of  each  convolution  function  to  be  zero, 
so  that  the  net  control  on  each  wall  is  exactly  zero  at  any  instant.  This  is  motivated  both 
by  physical  flow  control  devices  and  the  current  simulations  which  require  control  with 
zero  net  mass  flux.  We  will  seek  the  best  control  rule  to  interact  with  the  fluctuating 
part  of  the  turbulence  only. 

Note  also  that  the  weighting  functions  are  prescribed  at  the  outset  to  be  invariant 
in  time.  Though  the  method  used  requires  that  the  weights  be  optimized  by  considering 
finite  time  horizon  [0,T],  we  seek  to  approximate  the  steady^state  weights  at  the  “infinite 
time  horizon”  in  which  turbulent  fluctuations  near  the  wall  are  countered  by  a  fixed 
control  rule  at  the  wall  in  an  efficient  drag-reducing  manner. 

As  a  straightforward  extension  to  this  work,  one  may  also  optimize  nonlinear  control 
rules  of  a  form  similar  to  that  used  in  neural  networks,  which  may  be  written 


A 


A=1 


.  dui  du3 

CA  =  W^IA  *  - h  W2\  *  P  +  Wsx  *  {J,-—  +  Ha, 

on  on 

where  A  is  the  number  of  “nodes”  in  the  “hidden  layer”  of  the  network  and  g{^x)  is  an 


60 


T.  R.  Beviley  &  P.  Main 


Output  =  Control  Velocity 


Input  =  Flow  Measurements 

Figure  1.  Linear  Network.  The  flow  measurements  which  we  take  as  inputs  are  localized  mea¬ 
surements  on  the  wall  of  the  streamwise  drag,  the  pressure,  and  the  spanwise  drag.  The  flow 
measurements  are  convolved  with  the  weighting  functions  and  summed  to  determine  the  con¬ 
trol  The  input  flow  measurements  are  field  variables  and  are  indicated  with  heavy  lines — the 
corresponding  weights  are  convolution  functions  (in  the  continuous  case)  or  two  dimensional 
arrays  containing  a  stencil  of  weights  (in  the  discrete  case). 


Output  =  Control  Velocity 


Figure  2.  Nonlinear  Network.  The  output  of  several  simple  networks  A  similar  to  the  one 
depicted  in  Figure  1  (with  added  bias  weights  B\  connected  to  an  input  clamped  to  unity)  are 
used  as  the  arguments  to  activation  functions  g  at  the  hidden  nodes.  The  output  of  all  of 
the  hidden  nodes  ^(^a)  are  then  weighted  with  the  W\  and  summed  to  produce  the  control 


“activation  function”  which  will  be  prescribed.  Control  rules  of  this  form,  used  commonly 
in  neural  networks,  have  seen  a  broad  range  of  application  and  are  capable  of  representing 
very  general  nonlinear  relationships  (Hertz  et  al.^  1991).  The  task  at  hand  in  this  case  is 
to  optimize  the  weighting  functions  Wf^x  and  the  discrete  weights  Bx  and  Wx  such  that 
$  effectively  controls  the  flow  system.  This  type  of  network  is  illustrated  in  Figure  2. 

For  simplicity,  the  equations  necessary  to  optimize  the  linear  network  of  equation  (2.1) 
and  figure  1  will  be  derived  in  this  paper. 


3.  Governing  equations 

The  flow  system  we  consider  is  fully  developed  turbulent  channel  flow  with  periodic 
boundary  conditions  in  the  streamwise  (a^i)  and  spanwise  (2^3)  directions,  as  shown  in 
Figure  3.  Blowing  and  suction  through  the  computational  equivalent  of  holes  drilled 
in  the  wall  will  be  applied  according  to  the  linear  control  rule  (2.1).  The  control  rule 
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Figure  3.  Flow  configuration.  Blowing  and  suction  is  applied  through  closely  spaced  holes 
drilled  in  the  walls  to  control  the  flow.  This  control  will  be  coordinated  with  nearby  flow 
measurements  on  the  wall  in  a  manner  to  be  determined. 


optimized  for  this  configuration  should  also  be  effective  in  turbulent  boundary  layers  due 
to  the  similar  near- wall  behavior  of  these  flows. 

The  governing  equations  may  be  written  functionally  as 


Af{U)  =  0  in  (3.1a) 

with  boundary  conditions 

Ui  =  ^ni  ond^i^  (3.1b) 


where  ^  is  determined  by  the  linear  control  rule 


dui  duz 

=  wi  *  /i— - h  *  P  +  wz  * 

an  on 


(3.1c) 


and  prescribed  initial  conditions 


Ui  =  '«i(0) 


att  =  Q.  (3.1d) 


For  clarity,  all  differential  equations  will  be  written  in  operator  form  in  this  discussion, 
with  these  operators  defined  when  first  introduced.  The  (nonlinear)  Navier-Stokes  op¬ 
erator  for  the  present  case,  in  which  the  flow  is  assumed  to  have  uniform  density  and 
viscosity,  is  given  by 


/ 


M{U)  = 


4"  Uj 


dui 


dt  ’  ^  dxj 


d^ui  dp  1 

_j - - - 1 - 

p  OXi  p 


dxj 


duj 

dxj 


In  this  discussion,  Xi  is  the  stream  wise  direction,  X2  is  the  wall-normal  direction,  x^ 
is  the  spanwise  direction,  the  Ui^  are  the  corresponding  velocities,  p  is  the  pressure, 
p  is  the  density,  p  is  the  absolute  viscocity,  and  u  =  pijp  is  the  kinematic  viscocity. 
The  flow  is  sustained  by  a  mean  pressure  gradient  in  the  streamise  direction,  which  is 
modified  at  each  time  step  in  order  to  maintain  a  constant  mass  flux  through  the  channel. 
Define  also  n  as  a  wall-normal  unit  vector  directed  into  the  channel,  =  p  dui  f  dn\umii 
as  the  mean  skin  friction  on  the  wall  for  the  uncontrolled  channel  (averaged  in  space 
and  time),  Ur  ^  (^/p)^/^  as  the  mean  friction  velocity,  S  as  the  channel  half-width, 
and  Rcr  =  as  the  Reynolds  number  based  on  the  mean  friction  velocity  and  the 

channel  half  width.  The  flow  considered  in  this  work  is  taken  at  Rcr  =  180.  All  velocities 
are  normalized  by  the  friction  velocity  Ur^  and  therefore  may  also  be  marked  with  a  (*^) 
superscript.  All  lengths  are  normalized  by  S  unless  marked  with  a  (■*“)  superscript,  in 


f  K 
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which  case  they  are  normalized  by  the  wall  unit  vjur.  All  times  are  normalized  by  S/ur 
unless  marked  with  a  (+)  superscript,  in  which  case  they  are  normalized  by  vfu^. 

We  will  now  develop  a  systematic  procedure  to  optimize  the  weighting  functions 
The  first  step  is  to  define  the  control  problem  of  interest  mathematically  as  a  cost  func¬ 
tional  to  be  minimized. 


4.  Cost  functional 

The  objective  of  applying  the  control  in  this  problem  is  to  reduce  the  drag  without 
using  excessive  amounts  of  control  forcing.  Mathematically,  a  cost  functional  for  this 
problem  may  thus  be  expressed  as 

The  term  involving  /i  dui  / dn  is  the  drag  averaged  over  the  wall  F  and  the  time  interval 
of  interest  [0,T].  The  term  involving  is  an  expression  of  the  magnitude  of  the  control. 
These  two  terms  are  weighted  together  with  a  factor  P ,  which  represents  the  price  of  the 
control.  This  quantity  is  small  if  the  control  is  “cheap” ,  and  large  if  applying  the  control 
is  “expensive” . 

Minimization  of  J  corresponds  to  reducing  drag  while  maintaining  a  small  amount  of 
control  forcing. 


5.  Gradient  of  cost  functional 

We  now  develop  a  technique  to  compute  the  gradient  of  the  cost  functional  J  with 
respect  to  the  weighting  functions  w. 


5.1.  Perturbation  field 

Consider  first  the  Prechet  differential  of  the  flow  U  with  respect  to  w,  which  is  defined 
such  that 


•  1,.  U(w  ^  ew)  -  U{w) 

U  =  “7  lim - 

A  e-40  € 


--aL?. 


m{w) 


w^dT 


where  w  is  an  arbitrary  “update  direction”  to  the  weighting  function  w.  This  update 
direction  will  remain  undetermined  and  will  later  be  isolated  and  removed  from  the 
equation  for  the  perturbation  of  the  cost  functional  caused  by  a  perturbation  to  the 
control.  The  perturbation  field  U  is  governed  by  the  Frechet  differential  of  (3.1)  with 
respect  to  ty,  which  may  be  written: 


Alf  =  0 


in  Q. 


(5.1a) 


with  boundary  conditions 

Ui  =  ^Ui 


on 


(5.1b) 


where 
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•  dui  du3 

9  —  wi*  jjL— - W2  *  p  +  W3  *  [J'-r— 

on  an 


dui  dus 

A-Wi  *  fi— - h  +  , 

an  on 


and  with  initial  conditions 


Ui  =  0 


3.t  t  =  0. 


The  Frechet  differential  of  the  (non-linear)  Navier-Stokes  operator  is  given  by 


AU=^ 


diii  dui  ,  dui  d^Ui  1  dp 

dt  dxj  dxj  ^  dxj  p  dxi 


(5. Id) 


which  is  linear  in  the  perturbation  field  C/,  but  is  a  function  of  the  solution  U  of  the 
Navier-Stokes  problem,  so  that  A  =  A{U). 

The  Frechet  differential  of  the  cost  functional  J  with  respect  to  tt;  is: 

^  1,.  J{wAew)-J{w) 

J  =  —  lim - 

A  e-)-o  e 


Li: 


^J{w) 


Wk  dV 


It  is  seen  that  the  perturbation  of  the  cost  functional  is  a  function  of  the  perturbation 
of  the  flow  U .  The  linear  dependence  of  U  on  w  may,  in  theory,  be  found  directly  from 
(5.1).  However,  in  practice,  this  is  not  a  tractable  approach  due  to  the  excessively  large 
dimension  of  the  problem  under  consideration.  Thus,  we  seek  a  simpler  way  to  express 
the  above  equation  in  a  manner  in  which  the  gradient  S^J{w)ISiWf^  may  be  determined. 
It  is  for  this  reason  that  we  now  propose  the  definition  of  an  adjoint  field. 

5.2.  Definition  of  adjoint  field, 

As  discussed  in  previous  sections,  an  adjoint  operator  A*  may  be  defined  by  the  identity 

<AU,U>  =  <U,  A*U  >+b.  (5.3) 


Integration  by  parts  may  be  used  to  move  all  differential  operations  from  U  on  the  left 
hand  side  of  the  equation  to  U  on  the  right  hand  side,  resulting  in  an  expression  for  the 
adjoint  operator 


A*  U*  = 


(du*i  lap-’ 

^  \  dxj  dxi )  dx3  p  dxi 
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where  the  operator  A*  is  a  function  of  17,  and  an  expression  for  6,  which  contains  all  the 
boundary  terms: 


+ 


f  u'iuldV  -  [  u'i 

Jq  i=T  Ja 


■u*dV 


lt=o 


In  order  to  express  the  perturbation  of  the  cost  functional  in  a  usable  form,  we  now 
define  an  adjoint  state  by  the  system  of  equations 


A*U  =  0, 

with  mixed  boundary  conditions  on  the  walls 

Ui  f  ^  Wi  =  1  +  ^  ^  Wi 

—U2  n2  +  /  *  ^  ^W2 

Uz  +  /  *  w}3  =  $  *  103  , 

where 

du2 

f  =  p-2pu2  U2 

0X2 


(5.5a) 


(5.5b) 


and  with  initial  conditions 


b{x)  =  b{—x)^ 

Ui  =  0  at  t  =  T. 


(5.5c) 


5.3.  Identification  of  gradient 

Using  the  identity  (5.3)  and  the  definition  of  the  adjoint  in  (5.5),  we  can  algebraically 
manipulate  (5.3)  to  the  form 


A 


^J{w) 


Wk  dT 


A 


Gk,  dT , 


(5.6) 


where  is  some  function  of  the  solution  to  the  adjoint  problem  (5.5).  As  w  is  arbitrary, 
we  may  then  identify  the  expression  for  G^  as  ^  J{w) / It  may  be  shown  that  the 
resulting  expression  for  the  gradient  is 


^J(io) 

^J{w) 

^W2 

^J{w) 


=if 


dui 

*  /i— —  at 
on 


^pdt 


du3 

^  ^  at 

on 


Thus,  the  gradient  of  the  cost  ^7(io)  with  respect  to  modification  of  the  weights  w  of 
the  feedback  control  rule  has  been  represented  as  a  function  of  the  solution  of  an  adjoint 
problem.  This  adjoint  problem,  though  linear,  has  approximately  the  same  complexity 
as  the  Navier-Stokes  problem  itself,  as  can  be  seen  by  examining  Equation  (5.5)  and  the 
definition  of  A*.  Note  that  the  adjoint  field  evolves  backwards  over  the  domain  [0,T] 
from  initial  conditions  at  t  =  T]  this  is  the  natural  direction  in  time  to  march  these 
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equations  due  to  the  sign  of  the  viscous  and  time  derivatives  in  A*.  Note  also  that 
the  computation  of  the  adjoint  field  requires  the  storage  of  the  flow  field  over  the  entire 
interval  [0,T],  as  the  adjoint  operator  depends  on  the  flow  velocity,  i.e.  =  A*{U). 


6.  Gradient  update  to  control 

With  the  gradients  computed  using  the  adjoint  field,  a  control  rule  may  be  optimized 
using  a  gradient  algorithm  as  done  in  previous  sections. 


7.  Computational  Results 

At  the  time  this  document  went  to  press,  computational  results  were  not  yet  available. 
They  will  appear  in  the  version  of  to  be  submitted  to  J.  Fluid  Meek. 


8.  Discussion 

A  new  technique  for  optimizing  linear  and  nonlinear  feedback  control  rules  for  turbulent 
flows  has  been  presented.  This  technique  is  based  solely  on  the  equations  governing 
the  flow  and  a  mathematical  statement  of  the  control  objective,  thus  bypassing  the  ad 
hoc  identification  of  a  desired  state  trajectory  often  used  to  determine  feedback  control 
rules.  Also,  the  training  is  based  on  an  adjoint  (“sensitivity”)  field,  which  determines 
the  gradient  of  the  cost  with  respect  to  small  modifications  of  the  weights  in  a  rigorous 
manner.  Thus,  convergence  can  be  expected  to  be  much  better  than  for  an  reinforcement 
learning  approach  with  an  adaptive  algorithm. 

A  straightforward  extension  of  the  present  approach  is  to  take  into  account  past  mea¬ 
surements  in  the  control  rule.  Past  measurements,  which  may  easily  be  stored  in  an 
experimental  implementation,  may  give  additional  information  about  the  convection  ve¬ 
locity  of  flow  structures  which  cannot  be  determined  from  instantaneous  measurements 
alone.  It  is  also  possible  that  such  information  can  be  determined  by  recurrent  networks, 
in  which  the  inputs  of  the  control  network  include  the  outputs  of  the  network  from  the 
previous  time  step. 

Drawbacks  of  the  present  method  include  1)  an  accurate  mathematical  model  of  the 
flow  equations  and  boundary  conditions  are  needed  for  the  training,  and  2)  the  training 
algorithm  is  quite  complex,  requiring  simulation  on  a  supercomputer.  However,  this 
method  should  provide  insight  into  effective  new  control  rules  which  one  could  not  think 
of  otherwise,  and  which  can  be  further  modified  to  fit  practical  problems.  In  addition, 
they  may  be  used  to  guide  the  development  of  experimental  configurations,  revealing  the 
necessary  locations  of  sensors  with  respect  to  the  actuators  in  order  to  obtain  information 
relevant  to  effective  control  strategies. 
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